Juan Pablo Edwards Molina

Dataset

##      iso   ia    local  dose  cm
## 1 MES101 azox londrina 0.000 6.0
## 2 MES101 azox londrina 0.000 5.1
## 3 MES101 azox londrina 0.000 4.8
## 4 MES101 azox londrina 0.000 4.6
## 5 MES101 azox londrina 0.000 5.5
## 6 MES101 azox londrina 0.001 4.8
##        iso   ia   local dose  cm
## 235 MES202 tebu maringa    1 1.0
## 236 MES202 tebu maringa   10 0.1
## 237 MES202 tebu maringa   10 0.5
## 238 MES202 tebu maringa   10 0.8
## 239 MES202 tebu maringa   10 0.0
## 240 MES202 tebu maringa   10 0.0
iso1.ia1=subset(dat, iso==levels(iso)[1] & ia==levels(ia)[1])
dat1=within(iso1.ia1,
            {ctrl <- ave(cm, FUN=function(x) 1-(x/max(x)))}
            )
with(dat1, plot(ctrl~dose, pch=19, cex=0.7))

Ajustar o modelo

mod = glm(ctrl~ log(dose), data=subset(dat1, dose>0), family = quasibinomial (link="probit"))
summary(mod)
## 
## Call:
## glm(formula = ctrl ~ log(dose), family = quasibinomial(link = "probit"), 
##     data = subset(dat1, dose > 0))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.24271  -0.09083   0.03737   0.07073   0.17537  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.32541    0.03854   8.444 1.68e-08 ***
## log(dose)    0.17982    0.01014  17.739 6.49e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.01397955)
## 
##     Null deviance: 5.12238  on 24  degrees of freedom
## Residual deviance: 0.32573  on 23  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
par(mfrow=c(2,2)); plot(mod, which=1:3)
hnp(mod, xlab="Quantis teoricos", ylab="Residuos",main="Envelope GLM", cex.main=1, pch=4)
## Quasi-binomial model

layout(1)

Determinar o CE50

d50 = exp({log(0.5/(1-0.5)) - coef(mod)[[1]]} / coef(mod)[[2]])
d50
## [1] 0.1637196

Graficar o ajuste

plot(dat1$dose,dat1$ctrl, xlab = "dose",ylab = "prob",pch=19, cex=0.7)
ld <- seq(0, 10, 0.1)
lines(ld, predict(mod, data.frame(dose = ld), type = "response", col="red"))
abline(h=0.5, lty=3)
abline(v=d50[], lty=3, col="red")
text(d50[],0.001, "D50 = 0.163 ", adj = c(0,0), cex=0.8)