En un experimento descrito en Prentice (1976) se expuso una muestra de escarabajos a cierto pesticida. Tras cinco horas de exposición a distintos niveles de concentración del pesticida algunos de los escarabajos murieron y otros sobrevivieron. Los resultados para cada dosis aparecen en la tabla siguiente:

Dosis (\(\log_{10}CS_2 mgl^{-1}\))	N. insectos	N. muertos
1.6907	59	6
1.7242	60	13
1.7552	62	18
1.7842	56	28
1.8113	63	52
1.8369	59	53
1.8610	62	61
1.8839	60	60

Formula un modelo de regresión logística para analizar estos datos y estima la probabilidad de que muera un escarabajo expuesto durante cinco horas a una dosis de concentración 1.8.

Solución:

# Datos
dosis = c(rep(1.6907,59),rep(1.7242,60),rep(1.7552,62),rep(1.7842,56),
          rep(1.8113,63),rep(1.8369,59),rep(1.8610,62),rep(1.8839,60))
y = c(rep(1,6),rep(0,53),rep(1,13),rep(0,47),rep(1,18),rep(0,62-18),
      rep(1,28),rep(0,56-28),rep(1,52),rep(0,63-52),rep(1,53),rep(0,59-53),
      rep(1,61),rep(0,1),rep(1,60))

PropMuertes <- tapply(y,factor(dosis),mean)
plot(levels(factor(dosis)),PropMuertes,xlab="Dosis",ylab="Proporción de muertes")
lines(levels(factor(dosis)),PropMuertes)

# Ajuste del modelo
reg = glm(y~dosis,family='binomial')
summary(reg)

## 
## Call:
## glm(formula = y ~ dosis, family = "binomial")
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4922  -0.5986   0.2058   0.4512   2.3820  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -60.717      5.181  -11.72   <2e-16 ***
## dosis         34.270      2.912   11.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 645.44  on 480  degrees of freedom
## Residual deviance: 372.47  on 479  degrees of freedom
## AIC: 376.47
## 
## Number of Fisher Scoring iterations: 5

predict(reg,data.frame(dosis=1.8),type="response")

##         1 
## 0.7249464

dosisvec <- seq(1.69,1.89,0.001)
probMpred <- predict(reg,data.frame(dosis=dosisvec),type="response")
plot(levels(factor(dosis)),PropMuertes,xlab="Dosis",ylab="Proporción de muertes")
lines(levels(factor(dosis)),PropMuertes)
lines(dosisvec,probMpred,col="red")

Para aproximar el riesgo o tasa de clasificación errónea de la clasificación logística mediante leave-one-out:

library(boot)

## Warning: package 'boot' was built under R version 4.1.2

glm.fit <- glm(y~dosis,family='binomial',data=data.frame(cbind(dosis,y)))

cv.err <- cv.glm(data.frame(cbind(dosis,y)), glm.fit)
# Riesgo
cv.err$delta[1]

## [1] 0.1241303