Práctica 6: Regresión Logística

Author

Aritz Adin

Published

27/11/2024

1. Regresión Logística

Para estudiar la posible relación entre el número de satélites (cangrejos macho a su alrededor) y la anchura del caparazón de los cangrejos hembra tipo herradura, se realizó un experimento con 8 cangrejos hembra a los que se midió:

la anchura de los caparazones (\(x\)), medida en cm.
el número total de cangrejos macho que se utilizaron en el experimento para cada hembra (\(n\)).
el número de satélites o cangrejos macho que se situaron alrededor de cada hembra.

El objetivo es modelizar la proporción de satélites en función de la anchura del caparazón. Los datos se recogen en la siguiente tabla:

Anchura	Satélites	Total
23.25	5	14
23.75	4	14
24.75	17	28
25.75	21	39
26.75	15	22
27.75	20	24
28.75	15	18
29.25	14	14

Responde a las siguientes preguntas:

1. Introduce los datos en R mediante la función `data.frame`.

cangrejos <- data.frame(anchura=c(23.25,23.75,24.75,25.75,26.75,27.75,28.75,29.25),
                        satelites=c(5,4,17,21,15,20,15,14),
                        total=c(14,14,28,39,22,24,18,14))
print(cangrejos)

  anchura satelites total
1   23.25         5    14
2   23.75         4    14
3   24.75        17    28
4   25.75        21    39
5   26.75        15    22
6   27.75        20    24
7   28.75        15    18
8   29.25        14    14

2. Define una nueva variable llamada `prop` que indique la proporción de satélites para cada anchura.

cangrejos$prop <- cangrejos$satelites/cangrejos$total
print(cangrejos)

  anchura satelites total      prop
1   23.25         5    14 0.3571429
2   23.75         4    14 0.2857143
3   24.75        17    28 0.6071429
4   25.75        21    39 0.5384615
5   26.75        15    22 0.6818182
6   27.75        20    24 0.8333333
7   28.75        15    18 0.8333333
8   29.25        14    14 1.0000000

3. Ajusta el modelo probabilistico lineal \[\pi(x)=\alpha + \beta x\] donde \(y\) corresponde a la variable `prop` y \(x\) corresponde a la variable `anchura`

modelo1 <- lm(prop~anchura, weights=total, data=cangrejos)
summary(modelo1)


Call:
lm(formula = prop ~ anchura, data = cangrejos, weights = total)

Weighted Residuals:
     Min       1Q   Median       3Q      Max 
-0.39848 -0.28671 -0.00139  0.18904  0.60386 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.00206    0.39604  -5.055  0.00232 ** 
anchura      0.10081    0.01507   6.692  0.00054 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3641 on 6 degrees of freedom
Multiple R-squared:  0.8818,    Adjusted R-squared:  0.8621 
F-statistic: 44.78 on 1 and 6 DF,  p-value: 0.0005402

¿Cuáles son las estimaciones de los coeficientes? ¿Son estadísticamente significativos?

Las estimaciones de los coeficientes del modelo son: \(\hat{\alpha}=-2.002\) y \(\hat{\beta}=0.101\)

Ambos coeficientes son estadísticamente significativos, ya que los p-valores son pequeños.
¿Cuál es el coeficiente de determinación? Explica su significado.

\(R^2_\mbox{adjusted}=0.8621\). Es decir, el 86.21% de la variabilidad total es explicada por el modelo de regresión (variable`anchura).
¿El modelo ajustado es estadísticamente significativo?

Sí, el modelo es estadísticamente significativo. La función anova() presenta el análisis de la varianza con la descomposición de la variabilidad total en sus fuentes de variación y nos indica que la variable anchura es estadísticamente significativa.

4. ¿Cuáles son las predicciones de la proporción de satélites para cada valor de la anchura observada? ¿Y el número esperado de satélites por hembra?

data.frame(anchura=cangrejos$anchura,
           prop=cangrejos$prop,
           prop.pred=fitted(modelo1),
           satelite=cangrejos$satelites,
           satelite.pred=fitted(modelo1)*cangrejos$total)

  anchura      prop prop.pred satelite satelite.pred
1   23.25 0.3571429 0.3418064        5       4.78529
2   23.75 0.2857143 0.3922122        4       5.49097
3   24.75 0.6071429 0.4930236       17      13.80466
4   25.75 0.5384615 0.5938350       21      23.15957
5   26.75 0.6818182 0.6946465       15      15.28222
6   27.75 0.8333333 0.7954579       20      19.09099
7   28.75 0.8333333 0.8962694       15      16.13285
8   29.25 1.0000000 0.9466751       14      13.25345

5. Representa gráficamente las proporciones observadas y las predicciones lineales de las proporciones de satélites por cada anchura del caparazón.

plot(cangrejos$anchura, cangrejos$prop, xlab="Anchura (cm)", ylab="Prop", ylim=c(-0.1,1.1),
     main="Modelo probabilistico lineal")
abline(h=c(0,1), lty=2)
abline(modelo1, col="red", lwd=2)

6. Ajusta el modelo de regresión logística \[\mbox{logit}[\pi(x)] = \log \left(\frac{\pi(x)}{1-\pi(x)} \right) = \alpha + \beta x\] donde \(\pi\) corresponde a la variable `prop`, y \(x\) a la variable `anchura`.

modelo2 <- glm(prop~anchura, weights=total, family=binomial, data=cangrejos)
summary(modelo2)


Call:
glm(formula = prop ~ anchura, family = binomial, data = cangrejos, 
    weights = total)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -12.4421     2.7150  -4.583 4.59e-06 ***
anchura       0.5011     0.1051   4.767 1.87e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 34.0340  on 7  degrees of freedom
Residual deviance:  6.4542  on 6  degrees of freedom
AIC: 33.636

Number of Fisher Scoring iterations: 4

¿Cuáles son las estimaciones de los coeficientes? ¿Son estadísticamente significativos?

Las estimaciones de los coeficientes del modelo son: \(\hat{\alpha}=-12.442\) y \(\hat{\beta}=0.501\)

Ambos coeficientes son estadísticamente significativos, ya que los p-valores son pequeños.

Obtenemos la misma conclusión si interpretamos los intervalos de confianza de los parámetros del modelo:

confint(modelo2)

Waiting for profiling to be done...

                 2.5 %     97.5 %
(Intercept) -18.027915 -7.3310864
anchura       0.303846  0.7180261

¿El modelo ajustado es estadísticamente significativo?

7. ¿Cuál es la tasa incremental de cambio a los 26 cm? ¿Qué significa?

La tasa incremental de cambio en un punto \(x\) se calcula como \[\hat{\beta}*\hat{\pi}(x)*(1-\hat{\pi}(x))\]

Para una anchura de \(x=26\) cm su valor es:

beta <- coef(modelo2)[2]
hat.prop <- predict(modelo2, newdata=data.frame(anchura=26), type="response")

beta*hat.prop*(1-hat.prop)

 anchura 
0.115072

También se puede calcular mediante la función tasa.incremental()

source("tasa.incremental.R")
tasa.incremental(modelo2, 26, 27, decimales=4)

Argumentos de entrada:
  x1=26; x2=27; incremento=1 

Probabilidades estimadas:
  pi.x1=0.6427; pi.x2=0.7481; tasa incremental=0.1151 

Odds estimados:
 odds.x1=1.7991; odds.x2=2.9695; OR(odds.x2/odds.x1)=1.6506

Interpretación: Cuando la anchura del caparazón aumenta de 26 a 27 cm, la proporción de satélites por hembra aumenta en 0.1151.

8. ¿Cuál es la tasa incremental de cambio de probabilidad para pasar de 23 a 24 cm en la anchura del caparazón? ¿Y de 25 a 26 cm? ¿Qué significan?

tasa.incremental(modelo2, 23, 24, decimales=4)

Argumentos de entrada:
  x1=23; x2=24; incremento=1 

Probabilidades estimadas:
  pi.x1=0.2858; pi.x2=0.3977; tasa incremental=0.1023 

Odds estimados:
 odds.x1=0.4001; odds.x2=0.6604; OR(odds.x2/odds.x1)=1.6506

tasa.incremental(modelo2, 25, 26, decimales=4)

Argumentos de entrada:
  x1=25; x2=26; incremento=1 

Probabilidades estimadas:
  pi.x1=0.5215; pi.x2=0.6427; tasa incremental=0.1251 

Odds estimados:
 odds.x1=1.09; odds.x2=1.7991; OR(odds.x2/odds.x1)=1.6506

9. ¿Cuánto crece el odds estimado (odds ratio) por cada cm de aumento del caparazón? Proporciona un intervalo de confianza al 95%.

Cuando la anchura del caparazón aumenta en 1 cm, el odds ratio (\(OR=e^\beta\)) es de

exp(coef(modelo2)[2])

 anchura 
1.650584

Es decir, por cada incremento de 1 cm en la anchura del caparazón el odds estimado crece aproximadamente un 65%.

Para calcular un intervalo de confianza al 95% de \(e^\beta\), debemos exponenciar los límites inferior y superior del IC para el parámetro \(\beta\).

IC.beta <- confint(modelo2)[2,]
IC.beta

    2.5 %    97.5 % 
0.3038460 0.7180261

IC.OR <- exp(IC.beta)
IC.OR

   2.5 %   97.5 % 
1.355060 2.050382

Es decir, con un nivel de confianza del 95%, por cada incremento de 1 cm en la anchura del caparazón el odds estimado crece entre 1.36 y 2.05 veces (o lo que es lo mismo, crece entre un 35% y un 105%).

10. ¿Cuáles son las predicciones de la proporción de satélites para cada valor de la anchura observada?¿Y el número esperado de satélites por hembra?

data.pred <- data.frame(anchura=cangrejos$anchura,
                        prop=cangrejos$prop,
                        prop.pred=fitted(modelo2),
                        satelite=cangrejos$satelites,
                        satelite.pred=fitted(modelo2)*cangrejos$total)
data.pred

  anchura      prop prop.pred satelite satelite.pred
1   23.25 0.3571429 0.3119907        5      4.367870
2   23.75 0.2857143 0.3681261        4      5.153766
3   24.75 0.6071429 0.4902175       17     13.726090
4   25.75 0.5384615 0.6134872       21     23.926000
5   26.75 0.6818182 0.7237468       15     15.922429
6   27.75 0.8333333 0.8121823       20     19.492375
7   28.75 0.8333333 0.8771142       15     15.788056
8   29.25 1.0000000 0.9016724       14     12.623414

11. Representa gráficamente las proporciones observadas y las predicciones lineales y logísticas de las proporciones de satélites por cada anchura del caparazón.

plot(cangrejos$anchura, cangrejos$prop, xlab="Anchura (cm)", ylab="Prop", ylim=c(-0.1,1.1),
     main="Modelo regresión logística")
abline(h=c(0,1), lty=2)
abline(modelo1, col="red", lwd=2)
lines(data.pred$anchura, data.pred$prop.pred, col="blue", lwd=2)
legend(27, 0.3, legend=c("Modelo lineal","Regresión logística"), col=c("red", "blue"), lwd=2, bty="n")

12. Calcula la proporció estimada de cangrejos satélite para una hembra con anchura de caparazón de 28 cm. Calcula e interpreta un intervalo de confianza al 95% para la proporción verdadera.

data.pred <- data.frame(anchura=28)
data.pred <- cbind(data.pred,
                   predict(modelo2, data.pred, type="link", se.fit=TRUE))
data.pred

  anchura      fit   se.fit residual.scale
1      28 1.589535 0.290366              1

## Intervalo de confianza para el logit
alpha <- 0.05
lim.inf <- data.pred$fit-qnorm(1-alpha/2)*data.pred$se.fit
lim.sup <- data.pred$fit+qnorm(1-alpha/2)*data.pred$se.fit
IC.logit <- c(lim.inf,lim.sup)
IC.logit

[1] 1.020428 2.158642

## Intervalo de confianza para la probabilidad
IC.prob <- c(plogis(lim.inf),plogis(lim.sup))
IC.prob

[1] 0.7350560 0.8964736

2. Para entregar

Para estudiar cómo afecta una determinada droga sobre una especie de insectos, se realizó un experimento con 8 grupos de insectos a los que se les administró distintas dosis de dicha droga y se les midió:

la dosis administrada a cada grupo de insectos (\(x\))
el número total de insectos por grupo (\(n\))
el número de insectos muertos en cada grupo tras 5 horas de exposición a la droga

El objetivo es modelizar la proporción de insectos muertos en función de la dosis administrada. Los datos se recogen en la siguiente tabla:

Dosis	Insectos muertos	Total insectos
169.1	6	59
172.4	13	60
175.5	18	62
178.4	28	56
181.1	52	63
183.7	53	59
186.1	61	62
188.3	60	60

Responde a las siguientes preguntas:

2.1. Introduce los datos en R mediante la función `data.frame`.

2.2. Define una nueva variable llamada `prop` que indique la proporción de insectos muertos.

2.3. Ajusta el modelo probabilistico lineal \[\pi(x)=\alpha + \beta x\] donde \(y\) corresponde a la variable `prop` y \(x\) corresponde a la variable `dosis`

¿Cuáles son las estimaciones de los coeficientes? ¿Son estadísticamente significativos?
¿Cuál es el coeficiente de determinación? Explica su significado.
¿El modelo ajustado es estadísticamente significativo?

2.4. ¿Cuáles son las predicciones de la proporción de insectos muertos para cada dosis? ¿Y el número esperado de insectos muertos por dosis?

2.5. Representa gráficamente las proporciones observadas y las predicciones lineales de las proporciones de satélites por cada anchura del caparazón.

2.6. Ajusta el modelo de regresión logística \[\mbox{logit}[\pi(x)] = \log \left(\frac{\pi(x)}{1-\pi(x)} \right) = \alpha + \beta x\] donde \(\pi\) corresponde a la variable `prop`, y \(x\) a la variable `dosis`.

¿Cuáles son las estimaciones de los coeficientes? ¿Son estadísticamente significativos?
¿El modelo ajustado es estadísticamente significativo?

2.7. ¿Cuál es la tasa incremental de cambio para una dosis de 170? ¿Qué significa?

2.8. ¿Cuánto crece el odds estimado al aumentar la dosis del primer valor 169.1 al segundo valor 172.4?

2.9. ¿Cuáles son las predicciones de la proporción de satélites para cada valor de la anchura observada?

2.10. Representa gráficamente las proporciones observadas y las predicciones lineales y logísticas de las proporciones de satélites por cada anchura del caparazón.