У меня проблемы с построением кривой смертельной дозы с использованием модели биномиальной регрессии. Дело в том, что кривая ведет себя правильно, как показано ниже:
Идея будет выглядеть примерно так (рисунок 1).
Уведомлениемодели хорошо настроены (рисунок 2).
Проблема, с которой можно столкнуться (рисунок 3)
Ссылка envelr_bino.r: source ("http://www.ime.usp.br/~giapaula/envelr_bino")
Изображения
Рисунок 1: 
Рисунок 2: 
Рисунок 3: 
код:
resim = c(56,32,33,40,41)
resnao = c(4,28,27,20,19)
lnd = log(c(1,0.0265,0.053,0.0795,0.106)) # log of concentrations (put 1 on the first element for control when applying log 1 = 0)
dados = cbind(resim,resnao,lnd)
dados
dados = as.data.frame(dados)
attach(dados)
dados
dc = as.matrix(dados[,c(1,2)]) # as respostas
ajuste1 = glm(dc~lnd,family=binomial(link="logit"),data=dados)
ajuste1
anova(ajuste1,test="Chisq")
#ANOVA DO MODELO - temos que a
# a dose (lnd - logaritmo nepareriano da dose) é significativa
summary(ajuste1)
ajuste2 = glm(dc~lnd,family=binomial(link="probit"),data=dados)
ajuste2
anova(ajuste2,test="Chisq")
summary(ajuste2)
ajuste3 = glm(dc~lnd,family=binomial(link="cloglog"),data=dados)
ajuste3
anova(ajuste3,test="Chisq")
summary(ajuste3)
ajuste4 = glm(dc~lnd,family=binomial(link="cauchit"),data=dados)
ajuste4
anova(ajuste4,test="Chisq")
summary(ajuste4)
source("envelr_bino.R")
ntot = resim+resnao
x11()
par(mfrow=c(2,2))
envelr_bino(ajuste1)
envelr_bino(ajuste2)
envelr_bino(ajuste3)
envelr_bino(ajuste4)
x = seq(0,25,0.1)
m1 = exp(2.39+0.67*x)/(1+exp(2.39+0.67*x))
m2 = pnorm(1.42+0.039*x)
m3 = 1-exp(-exp(0.98+0.36*x))
m4 = 0.5+(atan(2.29+0.65*x)/pi)
par(mfrow=c(1,1))
x11()
plot(lnd,resim/(resim+resnao),pch=16,ylab="Mortes (%)",
xlab="ln(diluições)",xlim=c(-5,0),ylim=c(0,2),bty="n")
lines(x,m1, lty=2,lwd=2, col=2)
lines(x,m2, lty=3,lwd=2, col=4)
lines(x,m3, lty=6,lwd=2, col=3)
lines(x,m4, lty=1,lwd=2, col=1)
legend("top",lty=c(2,3,6,1),col=c(2,4,3,1),
lwd=2,c("logístico","probito","clog-log","cauchy"),bty="n")
lines(c(9,9),c(0,0.50),lty=3)
lines(c(0,9),c(0.50,0.50),lty=3)
legend(6,0.55,c("(9, 0.5)"),bty="n",cex=0.8)
РЕШЕНИЕ
renao = c(4,32,33,40,41)
resim = c(56,28,27,20,19)
lnd = c(0,0.0265,0.053,0.0795,0.106)
dados = cbind(resim,renao,lnd)
dados
dados = as.data.frame(dados)
dados
#cont <- ifelse(lnd==0, 1, 0)
#dados <- cbind(dados, cont)
dc = as.matrix(dados[,c(1,2)]) # as respostas
ajuste1 = glm(dc~lnd,family=binomial(link="logit"),data=dados)
ajuste1
anova(ajuste1,test="Chisq")
#ANOVA DO MODELO - temos que a
# a dose (lnd - logaritmo nepareriano da dose) é significativa
summary(ajuste1)
beta0 <- ajuste1$coefficients[1]
beta1 <- ajuste1$coefficients[2]
#model ajustado
xx <- seq(0, 0.15, 0.01)
eta <- beta0+beta1*xx
yy <- exp(eta)/(1+exp(eta))
x11()
plot(xx, yy, lwd=2, ylim=c(0,1), type='l')
points(lnd, resim/60, pch=16, col=2, cex=2)