Módulo 2. Unidad 2- Modelos Generalizados
1. Modelo de Regresión LogÃstico
Para el modelo de regresión logÃstico donde la respuesta es una
variable dicotómica se usará una base de datos de la librerÃa
mlbench con las variables: diabetes (positivo, negativo),
glucose (concentración de glucosa), pressure (presión diastólica mm Hg),
mass (IMC) y age (edad).
#install.packages("mlbench")
library("mlbench")
data("PimaIndiansDiabetes2", package = "mlbench")Para renombrar la base de datos:
data_diabetes<-PimaIndiansDiabetes2
names(data_diabetes)## [1] "pregnant" "glucose" "pressure" "triceps" "insulin" "mass" "pedigree"
## [8] "age" "diabetes"La estimación del modelo logÃstico teniendo como variable respuesta la diabetes (Si/No) se realiza a través de la función \(glm\), teniendo en cuenta la familia de la distribución binomial.
modelo2 <- glm( diabetes ~ glucose + pressure + mass + age, data = data_diabetes, family = binomial)
summary(modelo2)##
## Call:
## glm(formula = diabetes ~ glucose + pressure + mass + age, family = binomial,
## data = data_diabetes)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.673701 0.795070 -10.909 < 2e-16 ***
## glucose 0.034909 0.003523 9.910 < 2e-16 ***
## pressure -0.008317 0.008422 -0.988 0.323
## mass 0.092529 0.015315 6.042 1.52e-09 ***
## age 0.034352 0.008418 4.081 4.49e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 931.94 on 723 degrees of freedom
## Residual deviance: 695.17 on 719 degrees of freedom
## (44 observations deleted due to missingness)
## AIC: 705.17
##
## Number of Fisher Scoring iterations: 4Para obtener los OR´s y sus intervalos de confianza.
exp(coef(modelo2)) # exponenciar los coeficientes## (Intercept) glucose pressure mass age
## 0.0001710249 1.0355259404 0.9917173671 1.0969450393 1.0349486204exp(confint(modelo2)) # exponenciar los lÃmites del intervalo de confianza## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 3.415512e-05 0.0007744912
## glucose 1.028588e+00 1.0429135484
## pressure 9.753779e-01 1.0081878309
## mass 1.065070e+00 1.1310974425
## age 1.018122e+00 1.0523444507Para el análisis de los residuales de devianza, se puede obtener la tabla con la siguiente función:
anova(modelo2, test="Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: diabetes
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 723 931.94
## glucose 1 184.099 722 747.84 < 2.2e-16 ***
## pressure 1 3.224 721 744.61 0.07258 .
## mass 1 32.481 720 712.13 1.204e-08 ***
## age 1 16.964 719 695.17 3.809e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En la tabla anterior se observa como la devianza de los residuales del modelo nulo (sin covariables) es mayor respecto a los modelos con las covariables que se agregaron al modelo.
Para evaluar la bondad de ajuste del modelo se puede hacer uso de \(R^2\) de McFadden.
#install.packages("pscl")
library(pscl)
pR2(modelo2)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -347.5845435 -465.9686650 236.7682429 0.2540603 0.2789364 0.38529312. Modelo de Regresión Poisson
Para este modelo se utilizará una base de datos que se encuentra en la siguiente librerÃa:
#install.packages("ISwR")
library(ISwR)Los datos corresponden a casos de cáncer en diferentes ciudades, con las siguientes variables: ciudad, edad, población y casos.
Para cargar la base de datos:
data(eba1977)
datos_cancer = eba1977
head(datos_cancer)## city age pop cases
## 1 Fredericia 40-54 3059 11
## 2 Horsens 40-54 2879 13
## 3 Kolding 40-54 3142 4
## 4 Vejle 40-54 2520 5
## 5 Fredericia 55-59 800 11
## 6 Horsens 55-59 1083 6Para introducir la población como una variable offset en el modelo, se tranforma a logaritmo.
logpop = log(datos_cancer[ ,3]) #se crea la variable logaritmo de la población
datos_cancer2 = cbind(datos_cancer, logpop)
#se pega la nueva variable a la base de datos y se crea una nueva base.
datos_cancer2 ## city age pop cases logpop
## 1 Fredericia 40-54 3059 11 8.025843
## 2 Horsens 40-54 2879 13 7.965198
## 3 Kolding 40-54 3142 4 8.052615
## 4 Vejle 40-54 2520 5 7.832014
## 5 Fredericia 55-59 800 11 6.684612
## 6 Horsens 55-59 1083 6 6.987490
## 7 Kolding 55-59 1050 8 6.956545
## 8 Vejle 55-59 878 7 6.777647
## 9 Fredericia 60-64 710 11 6.565265
## 10 Horsens 60-64 923 15 6.827629
## 11 Kolding 60-64 895 7 6.796824
## 12 Vejle 60-64 839 10 6.732211
## 13 Fredericia 65-69 581 10 6.364751
## 14 Horsens 65-69 834 10 6.726233
## 15 Kolding 65-69 702 11 6.553933
## 16 Vejle 65-69 631 14 6.447306
## 17 Fredericia 70-74 509 11 6.232448
## 18 Horsens 70-74 634 12 6.452049
## 19 Kolding 70-74 535 9 6.282267
## 20 Vejle 70-74 539 8 6.289716
## 21 Fredericia 75+ 605 10 6.405228
## 22 Horsens 75+ 782 2 6.661855
## 23 Kolding 75+ 659 12 6.490724
## 24 Vejle 75+ 619 7 6.428105La estimación del modelo:
modelo3=glm(cases ~ city + age+ offset(logpop), family = poisson(link = "log"), data = datos_cancer2)
summary(modelo3)##
## Call:
## glm(formula = cases ~ city + age + offset(logpop), family = poisson(link = "log"),
## data = datos_cancer2)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.6321 0.2003 -28.125 < 2e-16 ***
## cityHorsens -0.3301 0.1815 -1.818 0.0690 .
## cityKolding -0.3715 0.1878 -1.978 0.0479 *
## cityVejle -0.2723 0.1879 -1.450 0.1472
## age55-59 1.1010 0.2483 4.434 9.23e-06 ***
## age60-64 1.5186 0.2316 6.556 5.53e-11 ***
## age65-69 1.7677 0.2294 7.704 1.31e-14 ***
## age70-74 1.8569 0.2353 7.891 3.00e-15 ***
## age75+ 1.4197 0.2503 5.672 1.41e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 129.908 on 23 degrees of freedom
## Residual deviance: 23.447 on 15 degrees of freedom
## AIC: 137.84
##
## Number of Fisher Scoring iterations: 5Para obtener los Riesgos relativos y sus intervalos de confianza:
exp(coef(modelo3)) # exponenciar los coeficientes## (Intercept) cityHorsens cityKolding cityVejle age55-59 age60-64
## 0.003581174 0.718880610 0.689667168 0.761612264 3.007213795 4.565884929
## age65-69 age70-74 age75+
## 5.857402508 6.403619032 4.135686847exp(confint(modelo3)) # exponenciar los lÃmites del intervalo de confianza## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 0.002373629 0.005212346
## cityHorsens 0.502694733 1.025912422
## cityKolding 0.475568043 0.995045687
## cityVejle 0.525131867 1.098950868
## age55-59 1.842951851 4.901008833
## age60-64 2.907180919 7.236296972
## age65-69 3.748295295 9.248885425
## age70-74 4.043044796 10.211923083
## age75+ 2.522891909 6.7624225723. Modelo de Regresión de Cox
Para la regresión de cox se utilizará una base de datos sobre un
ensayo clÃnico de dos tratamientos para cancer de pulmón en veteranos de
guerra de la librerÃa survival.
library(survival)
library(dplyr)
library(ggplot2)
library(ggfortify)Las variables contenidas en la base de datos son:
- trt: 1=standard 2=test
- celltype: 1=squamous, 2=small cell, 3=adeno, 4=large
- time: survival time in days
- status: censoring status
- karno: Karnofsky performance score (100=good)
- diagtime: months from diagnosis to randomization
- age: in years
- prior: prior therapy 0=no, 10=yes.
data(veteran)## Warning in data(veteran): data set 'veteran' not foundhead(veteran)## trt celltype time status karno diagtime age prior
## 1 1 squamous 72 1 60 7 69 0
## 2 1 squamous 411 1 70 5 64 10
## 3 1 squamous 228 1 60 3 38 0
## 4 1 squamous 126 1 60 9 63 10
## 5 1 squamous 118 1 70 11 65 10
## 6 1 squamous 10 1 20 5 49 0Inicialmente se analizará la tabla de sobrevida de Kaplan-Meier.
km_fit <- survfit(Surv(time, status) ~ 1, data=veteran)
#se debe incluir la variable tiempo y el evento/censura
summary(km_fit, times = c(1,30,60,90*(1:10)))## Call: survfit(formula = Surv(time, status) ~ 1, data = veteran)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 137 2 0.985 0.0102 0.96552 1.0000
## 30 97 39 0.700 0.0392 0.62774 0.7816
## 60 73 22 0.538 0.0427 0.46070 0.6288
## 90 62 10 0.464 0.0428 0.38731 0.5560
## 180 27 30 0.222 0.0369 0.16066 0.3079
## 270 16 9 0.144 0.0319 0.09338 0.2223
## 360 10 6 0.090 0.0265 0.05061 0.1602
## 450 5 5 0.045 0.0194 0.01931 0.1049
## 540 4 1 0.036 0.0175 0.01389 0.0934
## 630 2 2 0.018 0.0126 0.00459 0.0707
## 720 2 0 0.018 0.0126 0.00459 0.0707
## 810 2 0 0.018 0.0126 0.00459 0.0707
## 900 2 0 0.018 0.0126 0.00459 0.0707Gráfico de Kaplan-Meier:
autoplot(km_fit)
Gráfico de Kaplan-Meier diferenciado por tratamiento:
km_trt_fit <- survfit(Surv(time, status) ~ trt, data=veteran)
autoplot(km_trt_fit)
Para la estimación del modelo cox se utiliza la función \(coxph\):
cox <- coxph(Surv(time, status) ~ trt + celltype + karno + diagtime + age + prior , data = veteran)
summary(cox)## Call:
## coxph(formula = Surv(time, status) ~ trt + celltype + karno +
## diagtime + age + prior, data = veteran)
##
## n= 137, number of events= 128
##
## coef exp(coef) se(coef) z Pr(>|z|)
## trt 2.946e-01 1.343e+00 2.075e-01 1.419 0.15577
## celltypesmallcell 8.616e-01 2.367e+00 2.753e-01 3.130 0.00175 **
## celltypeadeno 1.196e+00 3.307e+00 3.009e-01 3.975 7.05e-05 ***
## celltypelarge 4.013e-01 1.494e+00 2.827e-01 1.420 0.15574
## karno -3.282e-02 9.677e-01 5.508e-03 -5.958 2.55e-09 ***
## diagtime 8.132e-05 1.000e+00 9.136e-03 0.009 0.99290
## age -8.706e-03 9.913e-01 9.300e-03 -0.936 0.34920
## prior 7.159e-03 1.007e+00 2.323e-02 0.308 0.75794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## trt 1.3426 0.7448 0.8939 2.0166
## celltypesmallcell 2.3669 0.4225 1.3799 4.0597
## celltypeadeno 3.3071 0.3024 1.8336 5.9647
## celltypelarge 1.4938 0.6695 0.8583 2.5996
## karno 0.9677 1.0334 0.9573 0.9782
## diagtime 1.0001 0.9999 0.9823 1.0182
## age 0.9913 1.0087 0.9734 1.0096
## prior 1.0072 0.9929 0.9624 1.0541
##
## Concordance= 0.736 (se = 0.021 )
## Likelihood ratio test= 62.1 on 8 df, p=2e-10
## Wald test = 62.37 on 8 df, p=2e-10
## Score (logrank) test = 66.74 on 8 df, p=2e-11cox_fit <- survfit(cox)
autoplot(cox_fit)
Para la validación del supuesto de proporcionalidad:
test.ph <- cox.zph(cox)
test.ph## chisq df p
## trt 0.2644 1 0.60712
## celltype 15.2274 3 0.00163
## karno 12.9352 1 0.00032
## diagtime 0.0129 1 0.90961
## age 1.8288 1 0.17627
## prior 2.1656 1 0.14113
## GLOBAL 34.5525 8 3.2e-05Por último, la validación a través de los de Schoenfeld se pueden obtener de la siguiente forma:
#install.packages("survminer")
library("survminer")
ggcoxzph(test.ph)