library(stringr)

source('read_mouse_liver_cancer_data.R')

NE = read_mouse_liver_cancer_data()

resp = cbind(NE$number_developed, NE$number_exposed - NE$number_developed)
months_on_study = NE$months_on_study
dose_amt = NE$dose_amt

## Fit a logit model with an interaction:
##
logit_model = glm(resp ~ months_on_study * dose_amt, family=binomial())
print(summary(logit_model))

## Fit a probit model with an interaction
##
probit_model = glm(resp ~ months_on_study * dose_amt, family=binomial(link=probit))
print(summary(probit_model))

## Fit a logit model without an interaction:
##
logit_model_2 = glm(resp ~ months_on_study + dose_amt, family=binomial())
print(summary(logit_model_2))

## Fit a probit model without an interaction
##
probit_model_2 = glm(resp ~ months_on_study * dose_amt, family=binomial(link=probit))
print(summary(probit_model_2))

if( FALSE ){
    ## Fit models withtout any factors:
    ## 
    months_on_study_number = as.numeric(str_replace(as.character(months_on_study), '[+]', '')) ## as a number

    ## Fit a logit model:
    ##
    logit_model = glm(resp ~ months_on_study_number * dose_amt, family=binomial())
    print(summary(logit_model))

    ## Fit a probit model:
    ##
    probit_model = glm(resp ~ months_on_study_number * dose_amt, family=binomial(link=probit))
    print(summary(probit_model))
}

## Estimate ED_01:
##
p = 0.01
logit_p = log(p/(1-p))
print(logit_p)

## Pick the largest month_of_study:
##
mask = months_on_study == '33+'
resp_33p = resp[mask, ]
dose_amt_33p = dose_amt[mask]

logit_model_3 = glm(resp_33p ~ dose_amt_33p, family=binomial())
print(summary(logit_model_3))

coefs = coefficients(logit_model_3)
ED_01 = ( logit_p - coefs[1] ) / coefs[2]
print(sprintf('ED_01 (all data)= %f', ED_01))

## Drop some of the low dose data:
##
resp_33p = resp_33p[-(1:2), ]
dose_amt_33p = dose_amt_33p[-(1:2)]

logit_model_4 = glm(resp_33p ~ dose_amt_33p, family=binomial())
print(summary(logit_model_4))

coefs = coefficients(logit_model_4)
ED_01 = ( logit_p - coefs[1] ) / coefs[2]
print(sprintf('ED_01 (only high dose data)= %f', ED_01))