Pregunta 2
A partir de la base drogadictos.csv:
library(readr)
datos <- read_csv("drogadictos.csv")
View(datos)
unique(datos$recaida)[1] 1 0datos$recaida <- factor(datos$recaida, levels = 0:1, labels = c("No","Sí"))
table(unclass(datos$recaida), datos$recaida)
No Sí
1 11 0
2 0 19datos$recaida = relevel(datos$recaida, ref=1)
unique(datos$terapia)[1] 0 1datos$terapia <- factor(datos$terapia, levels = 0:1, labels = c("Estándar","Nuevo"))
table(unclass(datos$terapia), datos$terapia)
Estándar Nuevo
1 15 0
2 0 15datos$terapia = relevel(datos$terapia, ref=1)
unique(datos$sexo)[1] 0 1datos$sexo <- factor(datos$sexo, levels = 0:1, labels = c("Masculino","Femenino"))
table(unclass(datos$sexo), datos$sexo)
Masculino Femenino
1 18 0
2 0 12table(datos$sexo,datos$recaida)
No Sí
Masculino 8 10
Femenino 3 9datos$sexo = relevel(datos$sexo, ref=1)
unique(datos$raza)[1] 1 0datos$raza <- factor(datos$raza, levels = 0:1, labels = c("Blanco","Afro"))
table(unclass(datos$raza), datos$raza)
Blanco Afro
1 18 0
2 0 12datos$terapia = relevel(datos$terapia, ref=1)
head(datos)# A tibble: 6 × 6
id tiempo recaida terapia sexo raza
<int> <int> <fctr> <fctr> <fctr> <fctr>
1 18 2 Sí Estándar Masculino Afro
2 13 2 Sí Estándar Masculino Afro
3 14 3 Sí Estándar Femenino Blanco
4 19 3 Sí Estándar Masculino Afro
5 25 4 Sí Estándar Femenino Afro
6 24 6 Sí Estándar Femenino AfroExploratorio
summary(datos) id tiempo recaida terapia sexo
Min. : 1.00 Min. : 2.00 No:11 Estándar:15 Masculino:18
1st Qu.: 8.25 1st Qu.: 7.00 Sí:19 Nuevo :15 Femenino :12
Median :15.50 Median : 9.50
Mean :15.50 Mean :11.73
3rd Qu.:22.75 3rd Qu.:16.00
Max. :30.00 Max. :24.00
raza
Blanco:18
Afro :12
library(GGally)
ggpairs(datos)
Planteamiento del modelo
Respuesta (D): Recaída Sí o No.
Exposición (E): Nueva terapia
Planteamos el modelo:
\[Logit[P(Recaida=Sí|Terapia)] =\\ \alpha + \beta*Terapia_{Nuevo}\]
lm1 <- glm(formula = recaida ~ terapia, family = binomial(link=log), data = datos)
summary(lm1)
Call:
glm(formula = recaida ~ terapia, family = binomial(link = log),
data = datos)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4823 -1.3537 0.9005 1.0108 1.0108
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4055 0.1826 -2.221 0.0264 *
terapiaNuevo -0.1054 0.2789 -0.378 0.7056
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 39.286 on 28 degrees of freedom
AIC: 43.286
Number of Fisher Scoring iterations: 5library(epitools)
epitab(datos$terapia, datos$recaida, method='riskratio')$tab
Outcome
Predictor No p0 Sí p1 riskratio lower upper p.value
Estándar 5 0.3333333 10 0.6666667 1.0 NA NA NA
Nuevo 6 0.4000000 9 0.6000000 0.9 0.5210192 1.554645 1
$measure
[1] "wald"
$conf.level
[1] 0.95
$pvalue
[1] "fisher.exact"exp(coef(lm1)) (Intercept) terapiaNuevo
0.6666667 0.9000000 exp(confint(lm1)) 2.5 % 97.5 %
(Intercept) 0.4154435 0.8652311
terapiaNuevo 0.4900046 1.5993025Prgeunta 3: Inclusión de raza y sexo
lm2 <- glm(data = datos, formula = recaida ~ terapia + sexo,
family = binomial(link=log), start = c(coef(lm1),0) )
lm3 <- glm(data = datos, formula = recaida ~ terapia + raza,
family = binomial(link=log), start = c(coef(lm1),0) )
summary(lm1); summary(lm2); summary(lm3)
Call:
glm(formula = recaida ~ terapia, family = binomial(link = log),
data = datos)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4823 -1.3537 0.9005 1.0108 1.0108
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.4055 0.1826 -2.221 0.0264 *
terapiaNuevo -0.1054 0.2789 -0.378 0.7056
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 39.286 on 28 degrees of freedom
AIC: 43.286
Number of Fisher Scoring iterations: 5
Call:
glm(formula = recaida ~ terapia + sexo, family = binomial(link = log),
data = datos, start = c(coef(lm1), 0))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6761 -1.2664 0.7505 1.0733 1.0909
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.57603 0.25547 -2.255 0.0241 *
terapiaNuevo -0.01904 0.26564 -0.072 0.9429
sexoFemenino 0.29442 0.27162 1.084 0.2784
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 38.222 on 27 degrees of freedom
AIC: 44.222
Number of Fisher Scoring iterations: 6
Call:
glm(formula = recaida ~ terapia + raza, family = binomial(link = log),
data = datos, start = c(coef(lm1), 0))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7829 -1.1855 0.6756 1.1246 1.3120
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.8606 0.2654 -3.243 0.00118 **
terapiaNuevo 0.2282 0.1613 1.415 0.15714
razaAfro 0.6324 0.2390 2.646 0.00814 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 34.697 on 27 degrees of freedom
AIC: 40.697
Number of Fisher Scoring iterations: 13lm4 <- glm(data = datos, formula = recaida ~ terapia + sexo + raza,
family = binomial(link=log) , start = c(coef(lm3),0))
summary(lm4)
Call:
glm(formula = recaida ~ terapia + sexo + raza, family = binomial(link = log),
data = datos, start = c(coef(lm3), 0))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9896 -1.1656 0.5453 0.9082 1.4730
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0849 0.3970 -2.733 0.00628 **
terapiaNuevo 0.4124 0.2421 1.703 0.08850 .
sexoFemenino 0.2637 0.2718 0.970 0.33199
razaAfro 0.6725 0.2528 2.660 0.00781 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 33.870 on 26 degrees of freedom
AIC: 41.87
Number of Fisher Scoring iterations: 23lm5 <- glm(data = datos, formula = recaida ~ terapia + sexo + raza + terapia:sexo,
family = binomial(link=log) , start = c(coef(lm4),0))
summary(lm5)
Call:
glm(formula = recaida ~ terapia + sexo + raza + terapia:sexo,
family = binomial(link = log), data = datos, start = c(coef(lm4),
0))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8457 -0.9999 0.6340 0.6681 1.4891
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.10866 0.39899 -2.779 0.00546 **
terapiaNuevo 0.22240 0.21682 1.026 0.30502
sexoFemenino 0.02141 0.31263 0.068 0.94540
razaAfro 0.88626 0.36057 2.458 0.01397 *
terapiaNuevo:sexoFemenino 0.64171 0.52461 1.223 0.22125
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 32.470 on 25 degrees of freedom
AIC: 42.47
Number of Fisher Scoring iterations: 6lm6 <- glm(data = datos, formula = recaida ~ terapia + sexo + raza + terapia:raza,
family = binomial(link=log) , start = c(coef(lm4),0))
summary(lm6)
Call:
glm(formula = recaida ~ terapia + sexo + raza + terapia:raza,
family = binomial(link = log), data = datos, start = c(coef(lm4),
0))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.9409 -1.0967 0.5744 0.9986 1.3158
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.86561 0.43693 -1.981 0.0476 *
terapiaNuevo 0.07143 0.46149 0.155 0.8770
sexoFemenino 0.36699 0.30928 1.187 0.2354
razaAfro 0.33367 0.42566 0.784 0.4331
terapiaNuevo:razaAfro 0.46052 0.54131 0.851 0.3949
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 33.284 on 25 degrees of freedom
AIC: 43.284
Number of Fisher Scoring iterations: 11exp(coef(lm1)) (Intercept) terapiaNuevo
0.6666667 0.9000000 exp(confint(lm1)) 2.5 % 97.5 %
(Intercept) 0.4154435 0.8652311
terapiaNuevo 0.4900046 1.5993025exp(coef(lm3)) (Intercept) terapiaNuevo razaAfro
0.4228969 1.2563569 1.8821425 exp(coef(lm4)) (Intercept) terapiaNuevo sexoFemenino razaAfro
0.3379398 1.5104120 1.3017355 1.9591388 Pregunta 3: Inclusión de raza y sexo en casos y controles
lm7 <- glm(data = datos, formula = recaida ~ terapia + raza + sexo + terapia:raza + terapia:sexo,
family = binomial(link=logit))
lm8 <- glm(data = datos, formula = recaida ~ terapia + raza + sexo,
family = binomial(link=logit))
library(lmtest)
lrtest(lm7, lm8)Likelihood ratio test
Model 1: recaida ~ terapia + raza + sexo + terapia:raza + terapia:sexo
Model 2: recaida ~ terapia + raza + sexo
#Df LogLik Df Chisq Pr(>Chisq)
1 6 -15.516
2 4 -16.685 -2 2.3388 0.3106summary(lm8)
Call:
glm(formula = recaida ~ terapia + raza + sexo, family = binomial(link = logit),
data = datos)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.2410 -1.0249 0.5690 0.8145 1.6141
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9855 0.9469 -1.041 0.2980
terapiaNuevo 0.6155 0.9375 0.657 0.5115
razaAfro 2.1088 1.0595 1.990 0.0465 *
sexoFemenino 1.3030 0.9175 1.420 0.1555
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 33.370 on 26 degrees of freedom
AIC: 41.37
Number of Fisher Scoring iterations: 4drop1(lm8, .~., test="LRT")Single term deletions
Model:
recaida ~ terapia + raza + sexo
Df Deviance AIC LRT Pr(>Chi)
<none> 33.370 41.370
terapia 1 33.816 39.816 0.4461 0.50420
raza 1 38.170 44.170 4.8000 0.02846 *
sexo 1 35.574 41.574 2.2040 0.13765
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lm9 <- glm(data = datos, formula = recaida ~ terapia + raza,
family = binomial(link=logit))
summary(lm9)
Call:
glm(formula = recaida ~ terapia + raza, family = binomial(link = logit),
data = datos)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.8542 -1.1919 0.6285 1.1231 1.2894
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.2594 0.7622 -0.340 0.7336
terapiaNuevo 0.3885 0.8922 0.435 0.6632
razaAfro 1.7809 0.9997 1.781 0.0748 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 35.574 on 27 degrees of freedom
AIC: 41.574
Number of Fisher Scoring iterations: 3drop1(lm9, .~., test="LRT")Single term deletions
Model:
recaida ~ terapia + raza
Df Deviance AIC LRT Pr(>Chi)
<none> 35.574 41.574
terapia 1 35.767 39.767 0.1928 0.66056
raza 1 39.286 43.286 3.7119 0.05403 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lm10 <- glm(data = datos, formula = recaida ~ terapia,
family = binomial(link=logit))
summary(lm10)
Call:
glm(formula = recaida ~ terapia, family = binomial(link = logit),
data = datos)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4823 -1.3537 0.9005 1.0108 1.0108
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6931 0.5477 1.266 0.206
terapiaNuevo -0.2877 0.7601 -0.378 0.705
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 39.429 on 29 degrees of freedom
Residual deviance: 39.286 on 28 degrees of freedom
AIC: 43.286
Number of Fisher Scoring iterations: 4exp(coef(lm8)) (Intercept) terapiaNuevo razaAfro sexoFemenino
0.3732669 1.8505853 8.2387288 3.6805011 exp(confint(lm8)) 2.5 % 97.5 %
(Intercept) 0.04729592 2.205997
terapiaNuevo 0.31337033 13.653478
razaAfro 1.23250078 89.291206
sexoFemenino 0.66959610 26.509498exp(coef(lm9)) (Intercept) terapiaNuevo razaAfro
0.7715257 1.4748184 5.9354567 exp(confint(lm9)) 2.5 % 97.5 %
(Intercept) 0.1518574 3.434028
terapiaNuevo 0.2658647 9.546525
razaAfro 0.9713270 55.740563exp(coef(lm10)) (Intercept) terapiaNuevo
2.00 0.75 exp(confint(lm10)) 2.5 % 97.5 %
(Intercept) 0.7107217 6.420991
terapiaNuevo 0.1627155 3.343129