Anova dos factores, modelos mixtos, regresión Logistica
Oscar Ramírez-Alán
6 de mayo de 2016
Como vimos anteriormente el anova de un factor trata de analizar el efecto de un bloque considerado mediante un modelo lineal. En este no se puede medir la interacción.
El Anova de dos vías es una técnica para el estudio de la relación entre una variable dependiente cuantitativa y dos variables cualitativas independientes. Por lo general, estamos interesados en saber si el nivel de la variable dependiente es diferente para los diferentes valores de las variables cualitativas.
ANOVA Factorial
Se trata de un diseño de bloques aleatorizados. Suponga se mide la producción, utilizando cuatro tipo de fertilizantes (cada finca es un bloque).
produc <- c(2.1, 2.2, 1.8, 2, 1.9, 2.2, 2.6, 2.7, 2.5, 2.8, 1.8, 1.9, 1.6, 2, 1.9, 2.1, 2, 2.2, 2.4, 2.1)
fert <- gl(4, 5)
finca <- factor(rep(1:5, 4))
#Permite crear una tabla de contingencia
xtabs(produc ~ finca + fert)## fert
## finca 1 2 3 4
## 1 2.1 2.2 1.8 2.1
## 2 2.2 2.6 1.9 2.0
## 3 1.8 2.7 1.6 2.2
## 4 2.0 2.5 2.0 2.4
## 5 1.9 2.8 1.9 2.1#permite crear un resumen para un factor
tapply(produc, fert, summary)## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.8 1.9 2.0 2.0 2.1 2.2
##
## $`2`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.20 2.50 2.60 2.56 2.70 2.80
##
## $`3`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.60 1.80 1.90 1.84 1.90 2.00
##
## $`4`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 2.10 2.10 2.16 2.20 2.40#permite crear un resumen para un factor
tapply(produc, finca, summary)## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.800 2.025 2.100 2.050 2.125 2.200
##
## $`2`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.900 1.975 2.100 2.175 2.300 2.600
##
## $`3`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.600 1.750 2.000 2.075 2.325 2.700
##
## $`4`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.000 2.000 2.200 2.225 2.425 2.500
##
## $`5`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.900 1.900 2.000 2.175 2.275 2.800#Recuerde tratar el ANOVA con los supuestos y estimaciones vistas en la clase anterior.aov(produc ~ finca + fert)->anofert
model.tables(anofert)#Permite ver el efecto para cada uno de los tratamientos## Tables of effects
##
## finca
## finca
## 1 2 3 4 5
## -0.090 0.035 -0.065 0.085 0.035
##
## fert
## fert
## 1 2 3 4
## -0.14 0.42 -0.30 0.02model.tables(anofert, type = "means")#Permite ver las medias de los tratamientos## Tables of means
## Grand mean
##
## 2.14
##
## finca
## finca
## 1 2 3 4 5
## 2.050 2.175 2.075 2.225 2.175
##
## fert
## fert
## 1 2 3 4
## 2.00 2.56 1.84 2.16Interacción
Se aplica con el objetivo de encontrar, si las variables dependientes, dependen de los factores (o niveles) que se están tratando.
#Su interpetación es: cuando existe interacción cuando el efecto de un factor depende del nivel del otro (o se comporte diferente para cada nivel del otro factor.#Recomiendo, siempre buscar la forma mas simple de interpretación
par(mfrow=c(1,2))
interaction.plot(fert, finca, produc, legend = T)
interaction.plot(finca, fert, produc, legend = T)
#permite crear un modelo de interacción
aov(produc ~ finca * fert)->anofert
#Observe que si hay un interacción
summary(anofert)## Df Sum Sq Mean Sq
## finca 4 0.088 0.0220
## fert 3 1.432 0.4773
## finca:fert 12 0.408 0.0340Otra forma de ver los resultados
#Creamos un modelo aditivo
lm(produc ~ finca + fert)->lmfert
#El modelo I y III (bajo la función "anova"), estima variación observada en las puntuaciones debidas al error experimental (Efectos Fijos).
anova(lmfert)## Analysis of Variance Table
##
## Response: produc
## Df Sum Sq Mean Sq F value Pr(>F)
## finca 4 0.088 0.02200 0.6471 0.6395716
## fert 3 1.432 0.47733 14.0392 0.0003137 ***
## Residuals 12 0.408 0.03400
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1(lmfert, test = "F") #tipo II ANOVA. Permite medir los efectos aleatorios (componentes de varianza).## Single term deletions
##
## Model:
## produc ~ finca + fert
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 0.408 -61.844
## finca 4 0.088 0.496 -65.938 0.6471 0.6395716
## fert 3 1.432 1.840 -37.719 14.0392 0.0003137 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diseño de bloques balanceado incompleto
#Utilizaremos la base de datos rabbit
#?rabbit, utilice este comando para revisar la información de la base de datos
library(faraway)
head(rabbit)## treat gain block
## 1 f 42.2 b1
## 2 b 32.6 b1
## 3 c 35.2 b1
## 4 c 40.9 b2
## 5 a 40.1 b2
## 6 b 38.1 b2#Visualizamos el modelo no balanceado
matrix(rabbit$treat,nrow=3)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] "f" "c" "c" "a" "e" "b" "d" "a" "d" "f"
## [2,] "b" "a" "f" "e" "c" "f" "a" "e" "b" "d"
## [3,] "c" "b" "d" "c" "d" "e" "b" "f" "e" "a"matrix(rabbit$block,nrow=3)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] "b1" "b2" "b3" "b4" "b5" "b6" "b7" "b8" "b9" "b10"
## [2,] "b1" "b2" "b3" "b4" "b5" "b6" "b7" "b8" "b9" "b10"
## [3,] "b1" "b2" "b3" "b4" "b5" "b6" "b7" "b8" "b9" "b10"#nótese que los bloques tienen misma cantidad de repeticiones, no así los tratamientos
xtabs(gain~treat+block, data=rabbit) #Note los faltantes de datos que son identificados con "0.0".## block
## treat b1 b10 b2 b3 b4 b5 b6 b7 b8 b9
## a 0.0 37.3 40.1 0.0 44.9 0.0 0.0 45.2 44.0 0.0
## b 32.6 0.0 38.1 0.0 0.0 0.0 37.3 40.6 0.0 30.6
## c 35.2 0.0 40.9 34.6 43.9 40.9 0.0 0.0 0.0 0.0
## d 0.0 42.3 0.0 37.5 0.0 37.3 0.0 37.9 0.0 27.5
## e 0.0 0.0 0.0 0.0 40.8 32.0 40.5 0.0 38.5 20.6
## f 42.2 41.7 0.0 34.3 0.0 0.0 42.8 0.0 51.9 0.0#Exploramos un poco los datos
plot(gain~ block + treat, rabbit, col=456)
#Que podemos concluir?
#Creamos un modelo lineal de regresión
lm.rabbit <- lm(gain~ block+treat, data=rabbit)
anova(lm.rabbit)## Analysis of Variance Table
##
## Response: gain
## Df Sum Sq Mean Sq F value Pr(>F)
## block 9 730.39 81.154 8.0738 0.0002454 ***
## treat 5 158.73 31.745 3.1583 0.0381655 *
## Residuals 15 150.77 10.052
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Observamos que tanto el bloque como el tratamiento son significativos.#Hasta aquí, se hace una diferencia, porque el diseño no es ortogonal a causa del carácter incompleto. Que tabla es adecuado para probar el efecto del tratamiento o el efecto de bloque?
summary(lm.rabbit)##
## Call:
## lm(formula = gain ~ block + treat, data = rabbit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.9583 -1.6146 -0.6083 1.9396 4.3028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.01389 2.58863 13.912 5.59e-10 ***
## blockb10 3.29722 2.79604 1.179 0.25667
## blockb2 4.13333 2.69433 1.534 0.14583
## blockb3 -1.80278 2.69433 -0.669 0.51360
## blockb4 8.79444 2.79604 3.145 0.00667 **
## blockb5 2.30556 2.79604 0.825 0.42253
## blockb6 5.40833 2.69433 2.007 0.06309 .
## blockb7 5.77778 2.79604 2.066 0.05651 .
## blockb8 9.42778 2.79604 3.372 0.00419 **
## blockb9 -7.48056 2.79604 -2.675 0.01729 *
## treatb -1.74167 2.24182 -0.777 0.44930
## treatc 0.40000 2.24182 0.178 0.86078
## treatd 0.06667 2.24182 0.030 0.97667
## treate -5.22500 2.24182 -2.331 0.03413 *
## treatf 3.30000 2.24182 1.472 0.16169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.17 on 15 degrees of freedom
## Multiple R-squared: 0.855, Adjusted R-squared: 0.7197
## F-statistic: 6.318 on 14 and 15 DF, p-value: 0.000518#Probamos solo analizar el bloque incompleto, para comparar el modelo
lm.rabbit2<- lm(gain ~treat,data=rabbit)
summary(lm.rabbit2)##
## Call:
## lm(formula = gain ~ treat, data = rabbit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.880 -3.050 1.200 2.825 9.320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.300 2.494 16.960 7.28e-15 ***
## treatb -6.460 3.527 -1.831 0.0795 .
## treatc -3.200 3.527 -0.907 0.3733
## treatd -5.800 3.527 -1.644 0.1131
## treate -7.820 3.527 -2.217 0.0363 *
## treatf 0.280 3.527 0.079 0.9374
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.577 on 24 degrees of freedom
## Multiple R-squared: 0.2821, Adjusted R-squared: 0.1326
## F-statistic: 1.886 on 5 and 24 DF, p-value: 0.1342coef(summary(lm.rabbit2))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.30 2.494173 16.95952784 7.282300e-15
## treatb -6.46 3.527294 -1.83143247 7.947811e-02
## treatc -3.20 3.527294 -0.90721113 3.733145e-01
## treatd -5.80 3.527294 -1.64432018 1.131471e-01
## treate -7.82 3.527294 -2.21699720 3.634801e-02
## treatf 0.28 3.527294 0.07938097 9.373878e-01anova(lm.rabbit, lm.rabbit2)## Analysis of Variance Table
##
## Model 1: gain ~ block + treat
## Model 2: gain ~ treat
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 15 150.77
## 2 24 746.51 -9 -595.74 6.5854 0.0007602 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Concluimos el análisis de bloques es bien para este caso. Notar que hay significancia entre los modelos.#Otra manera de probar, que recomiendo, cuando el modelo es inbalanceado es utilizar la función "drop1", del modelo que se selecciona.
drop1(lm.rabbit, test="F") ## Single term deletions
##
## Model:
## gain ~ block + treat
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 150.77 78.437
## block 9 595.74 746.51 108.426 6.5854 0.0007602 ***
## treat 5 158.73 309.50 90.013 3.1583 0.0381655 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Asume que todos los efectos están contenido en el modeloAjuste de un modelo lineal mixto (bloques).
Aunque el modelo lineal no es necesariamente inadecuado, la inclusión del bloque como un factor fijo puede también utilizarse. En este sentido el modelo mixto equivalente seria mucho más parsimonioso en término del número de parámetros
# nos interesa analizar el efecto del factor (dieta) sobre la respuesta (ganancia de peso) una vez que se haya tenido en cuenta la variabilidad atribuible al bloque (camada).
library(nlme)
lme.rabbit1 <- lme(gain~ treat, random=~1|block, data=rabbit)
anova(lme.rabbit1)## numDF denDF F-value p-value
## (Intercept) 1 15 590.5297 <.0001
## treat 5 15 3.2818 0.0336summary(lme.rabbit1)## Linear mixed-effects model fit by REML
## Data: rabbit
## AIC BIC logLik
## 166.3569 175.7813 -75.17844
##
## Random effects:
## Formula: ~1 | block
## (Intercept) Residual
## StdDev: 4.657853 3.175519
##
## Fixed effects: gain ~ treat
## Value Std.Error DF t-value p-value
## (Intercept) 39.53540 2.130338 15 18.558278 0.0000
## treatb -2.50718 2.208700 15 -1.135139 0.2741
## treatc -0.18407 2.208700 15 -0.083340 0.9347
## treatd -0.88516 2.208700 15 -0.400759 0.6942
## treate -5.64602 2.208700 15 -2.556263 0.0219
## treatf 2.81003 2.208700 15 1.272254 0.2227
## Correlation:
## (Intr) treatb treatc treatd treate
## treatb -0.518
## treatc -0.518 0.500
## treatd -0.518 0.500 0.500
## treate -0.518 0.500 0.500 0.500
## treatf -0.518 0.500 0.500 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.3794203 -0.5213453 -0.1701517 0.6456859 1.4148990
##
## Number of Observations: 30
## Number of Groups: 10Otros métodos aplicables, a los modelos mixtos
Comparamos modelos utilizando la función REML (Restricted / residual / reduced maximum likelihood). a) modelo sin efectos aleatorios (utilizando la función “gls”) b) modelo que incluya la varianza aleatoria de la constante (gran media), dentro de cada bloque (utilizando la función “lme”)
#Comparamos los modelos.
#un modelo sin efecto aleatorio
lme.rabbit0 <- gls(gain~ treat, data=rabbit) #Se utilza gls para poder compararlo, el resultados es similar a un modelo lineal
#un modelo con efecto aleatorio
lme.rabbit1 <- lme(gain~ treat, data=rabbit, random=~1|block)
#En ambos fue utilizando un estimador REML (bajo la función gls, que permite comparar los parámetros de modelos. Si se hubiera hecho con lm, no es posible compararlos)
anova(lme.rabbit0, lme.rabbit1)## Model df AIC BIC logLik Test L.Ratio p-value
## lme.rabbit0 1 7 174.2621 182.5085 -80.13107
## lme.rabbit1 2 8 166.3569 175.7813 -75.17844 1 vs 2 9.905255 0.0016#Vemos que el modelo con el efecto aleatorio bloque es más parsimonioso que el modelo sin el efecto aleatorio.#mejor opcion comparar modelos anidados como en el caso de los efectos
#aleatorios.
lme.rabbit2 <- lme(gain~ 1, data=rabbit, random=~1|block, method="ML")
lme.rabbit3 <- lme(gain~ treat, data=rabbit, random=~1|block, method="ML")
anova(lme.rabbit2, lme.rabbit3)## Model df AIC BIC logLik Test L.Ratio p-value
## lme.rabbit2 1 3 188.8307 193.0343 -91.41536
## lme.rabbit3 2 8 183.7730 194.9825 -83.88648 1 vs 2 15.05776 0.0101#Y, de nuevo, vemos que el modelo con el efecto fijo del tratamiento es más parsimonioso que el modelo que contiene únicamente el efecto de la constante.Analice los siguientes datos
Determine si el nivel enzimático depende del genotipo y el sexo.
enzyma<-read.table("Enzimas.txt", header=T)
head(enzyma,15)## Enzyma Sex Genotype
## 1 2.838 Female FF
## 2 4.216 Female FF
## 3 2.889 Female FF
## 4 4.198 Female FF
## 5 3.550 Female FS
## 6 4.556 Female FS
## 7 3.087 Female FS
## 8 1.943 Female FS
## 9 3.620 Female SS
## 10 3.079 Female SS
## 11 3.586 Female SS
## 12 2.669 Female SS
## 13 1.884 Male FF
## 14 2.283 Male FF
## 15 4.939 Male FFaov(Enzyma~Sex+Genotype, data=enzyma)->ano1
summary(ano1)## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.357 0.3567 0.523 0.478
## Genotype 2 0.484 0.2419 0.355 0.706
## Residuals 20 13.638 0.6819aov(Enzyma~Sex*Genotype, data=enzyma)->ano2
summary(ano2)## Df Sum Sq Mean Sq F value Pr(>F)
## Sex 1 0.357 0.3567 0.489 0.493
## Genotype 2 0.484 0.2419 0.332 0.722
## Sex:Genotype 2 0.512 0.2558 0.351 0.709
## Residuals 18 13.126 0.7292aov(Enzyma~Genotype*Sex, data=enzyma)->ano3
summary(ano3)## Df Sum Sq Mean Sq F value Pr(>F)
## Genotype 2 0.484 0.2419 0.332 0.722
## Sex 1 0.357 0.3567 0.489 0.493
## Genotype:Sex 2 0.512 0.2558 0.351 0.709
## Residuals 18 13.126 0.7292model.tables(ano1, type="means")## Tables of means
## Grand mean
##
## 3.230667
##
## Sex
## Sex
## Female Male
## 3.353 3.109
##
## Genotype
## Genotype
## FF FS SS
## 3.342 3.030 3.320model.tables(ano1, type="effects")## Tables of effects
##
## Sex
## Sex
## Female Male
## 0.12192 -0.12192
##
## Genotype
## Genotype
## FF FS SS
## 0.11096 -0.20042 0.08946#Comparamos modelos anteriores
lme.enzima <- lme(Enzyma~Sex, random=~1|Genotype, data=enzyma)
summary(lme.enzima)## Linear mixed-effects model fit by REML
## Data: enzyma
## AIC BIC logLik
## 65.64999 70.01415 -28.82499
##
## Random effects:
## Formula: ~1 | Genotype
## (Intercept) Residual
## StdDev: 1.930523e-05 0.8011862
##
## Fixed effects: Enzyma ~ Sex
## Value Std.Error DF t-value p-value
## (Intercept) 3.352583 0.2312825 20 14.495618 0.0000
## SexMale -0.243833 0.3270829 20 -0.745479 0.4647
## Correlation:
## (Intr)
## SexMale -0.707
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.75937046 -0.59453512 -0.09766768 0.41001705 2.28442527
##
## Number of Observations: 24
## Number of Groups: 3lme.enzima2 <- lme(Enzyma~Genotype, random=~1|Sex, data=enzyma)
summary(lme.enzima2)## Linear mixed-effects model fit by REML
## Data: enzyma
## AIC BIC logLik
## 67.31097 72.53358 -28.65548
##
## Random effects:
## Formula: ~1 | Sex
## (Intercept) Residual
## StdDev: 2.552969e-05 0.8163409
##
## Fixed effects: Enzyma ~ Genotype
## Value Std.Error DF t-value p-value
## (Intercept) 3.341625 0.2886201 20 11.577936 0.0000
## GenotypeFS -0.311375 0.4081705 20 -0.762855 0.4545
## GenotypeSS -0.021500 0.4081705 20 -0.052674 0.9585
## Correlation:
## (Intr) GntyFS
## GenotypeFS -0.707
## GenotypeSS -0.707 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.78555916 -0.62167655 -0.01071856 0.43467593 1.95674990
##
## Number of Observations: 24
## Number of Groups: 2lme.enzima3 <- lme(Enzyma~ 1, data=enzyma, random=~1|Sex, method="ML")
lme.enzima4 <- lme(Enzyma~ Genotype, data=enzyma, random=~1|Sex, method="ML")
anova(lme.enzima3, lme.enzima4)## Model df AIC BIC logLik Test L.Ratio p-value
## lme.enzima3 1 3 61.97973 65.51390 -27.98987
## lme.enzima4 2 5 65.16398 71.05425 -27.58199 1 vs 2 0.8157533 0.6651interaction.plot(enzyma$Genotype,enzyma$Sex,enzyma$Enzyma)
#Se utiliza la técnica del Análisis de Varianza Anidado (ANOVA anidada) cuando se tiene una variable de medición y dos o más variables nominales. Las variables nominales se anidan, lo que significa que cada valor de una variable nominal (los subgrupos) se encuentra en combinación con un solo valor de la variable nominal de más alto nivel (los grupos). El nivel nominal de la variable superior puede ser Modelo I o modelo II pero el nivel nominal de las variables más bajos deben ser el Modelo II.
#El análisis de varianza anidado es una extensión del análisis de varianza de una vía en que se divide cada grupo en subgrupos.
# En teoría, estos subgrupos se eligen al azar de un conjunto más amplio posible de los subgrupos. Un análisis de varianza anidado tiene una hipótesis nula para cada nivel. En el ANOVA anidada, una hipótesis nula sería que los subgrupos dentro de cada grupo tienen promedios iguales, la segunda hipótesis nula sería que los grupos tienen los mismas medias.Regresión logística
# Son una extensión de los modelos lineales que permiten utilizar distribuciones no normales de los errores (binomiales, Poisson, Gamma, etc), y varianzas no constantes.
#Si tenemos una variable que describe una respuesta en forma de dos posibles eventos (0 y 1) y queremos estudiar el efecto que otras variables independientes tienen sobre ella, el modelo de regresión logística binaria puede resultar de gran utilidad.
#El objetivo primordial que resuelve esta técnica es el de modelar cómo influye en la probabilidad de aparición de un suceso, habitualmente dicotómico, la presencia o no de diversos factores y el valor o nivel de los mismos.
#De todos es sabido que este tipo de situaciones se aborda mediante técnicas de regresión. Sin embargo, la metodología de la regresión lineal no es aplicable ya que ahora la variable respuesta sólo presenta dos valores (nos centraremos en el caso dicotómico), como puede ser presencia/ausencia.
#La regresión logística se utiliza para explicar y predecir una variable categórica binaria en función de las variables independientes (cuantitativas o cualitativas).
#algunos supuestos: a) La función de la regresión es no lineal, b) sigue un comportamiento binomial, por lo que invalida el supuesto de normalidad c) la varianza de la variable dicotómica no es constante.
#Los GLM permiten permiten especificar distintos tipos de distribución de errores:
#Poisson, muy útiles para conteos de acontecimientos .Ejemplo: número de heridos por accidentes de tráfico; número de hogares asegurados que dan parte de siniestro al día.
# Binomiales, de gran utilidad para proporciones y datos de presencia/ausencia Ejemplo: tasas de mortalidad; tasas de infección; porcentaje de siniestros mortales.
# Gamma, muy útiles con datos que muestran un coeficiente de variación constante, esto es, en donde la varianza aumenta según aumenta la media de la muestra de manera constante Ejemplo: número de heridos en función del número de siniestros
# Exponenciales, muy útiles para los análisis de supervivencia #Los datos "plasma" correspondiente al paquete "HSAUR" contiene (29) observaciones de tres variables: "fibrinogen" (Fibrinógeno), "globulin" ( Globulina) y "ESR" (Erythrocyte sedimentation rate).
data("plasma", package = "HSAUR")
head(plasma, 8)## fibrinogen globulin ESR
## 1 2.52 38 ESR < 20
## 2 2.56 31 ESR < 20
## 3 2.19 33 ESR < 20
## 4 2.18 31 ESR < 20
## 5 3.41 37 ESR < 20
## 6 2.46 36 ESR < 20
## 7 3.22 38 ESR < 20
## 8 2.21 37 ESR < 20par(mfrow=c(1,2))
#generamos gráficos de densidades para observar la distribución condicional de una variable, respecto a los cambios categóricos. "Solo funciona con variables categóricas".
cdplot(ESR ~ fibrinogen, data = plasma)
cdplot(ESR ~ globulin, data = plasma)
#Generamos el modelo, utilizando una función logística, dado el parámetro "binomial".
plasma_glm_1 <- glm(ESR ~ fibrinogen , data = plasma, family = binomial())
summary(plasma_glm_1)##
## Call:
## glm(formula = ESR ~ fibrinogen, family = binomial(), data = plasma)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9298 -0.5399 -0.4382 -0.3356 2.4794
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.8451 2.7703 -2.471 0.0135 *
## fibrinogen 1.8271 0.9009 2.028 0.0425 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 30.885 on 31 degrees of freedom
## Residual deviance: 24.840 on 30 degrees of freedom
## AIC: 28.84
##
## Number of Fisher Scoring iterations: 5#Con esto podemos concluir que la variable fibrinógeno es apropiada para explicar ESR.#Ajustamos un segundo modelo incluyendo ambas variables exploratorias.
plasma_glm_2 <- glm(ESR ~ fibrinogen + globulin, data = plasma, family = binomial())
summary(plasma_glm_2)##
## Call:
## glm(formula = ESR ~ fibrinogen + globulin, family = binomial(),
## data = plasma)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9683 -0.6122 -0.3458 -0.2116 2.2636
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.7921 5.7963 -2.207 0.0273 *
## fibrinogen 1.9104 0.9710 1.967 0.0491 *
## globulin 0.1558 0.1195 1.303 0.1925
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 30.885 on 31 degrees of freedom
## Residual deviance: 22.971 on 29 degrees of freedom
## AIC: 28.971
##
## Number of Fisher Scoring iterations: 5#Note que la variable "globulin", es no significativa, por lo que concluimos que no está asociada con los niveles de ESR.#Comparamos los dos modelos usando la función "anova".
anova(plasma_glm_1, plasma_glm_2, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: ESR ~ fibrinogen
## Model 2: ESR ~ fibrinogen + globulin
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 30 24.840
## 2 29 22.971 1 1.8692 0.1716library(ggplot2)
ggplot(plasma, aes(x=fibrinogen, y=ESR)) + geom_point()+ stat_smooth(method="glm")
ggplot(plasma, aes(x=globulin, y=ESR)) + geom_point()+ stat_smooth(method="glm")
Analice los siguientes datos, para determinar si la sobrevivencia depende de los niveles de hemoglobina y bilirubina en 70 individuos
enfermedad<-read.csv("enfermedad.csv", header=T, sep=";")
head(enfermedad)## sobrevivencia hemo bili
## 1 1 18.7 2.2
## 2 1 17.8 2.7
## 3 1 17.8 2.5
## 4 1 17.6 4.1
## 5 1 17.6 3.2
## 6 1 17.6 1.0enfermedad_glm_1<-glm(sobrevivencia~hemo+bili, family = binomial(),data=enfermedad)
summary(enfermedad_glm_1)##
## Call:
## glm(formula = sobrevivencia ~ hemo + bili, family = binomial(),
## data = enfermedad)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.93272 0.02706 0.09990 0.32891 2.30324
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8086 2.4853 -0.325 0.74490
## hemo 0.5616 0.1730 3.246 0.00117 **
## bili -0.7851 0.4212 -1.864 0.06236 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.056 on 69 degrees of freedom
## Residual deviance: 31.090 on 67 degrees of freedom
## AIC: 37.09
##
## Number of Fisher Scoring iterations: 7#Podemos concluir que la hemoglobila es una factor pronóstico importante.
cdplot(as.factor(sobrevivencia) ~ hemo, data = enfermedad)
cdplot(as.factor(sobrevivencia) ~ bili, data = enfermedad)
#Escribir la funcion
coef(enfermedad_glm_1)## (Intercept) hemo bili
## -0.8086451 0.5615559 -0.7850649#Estima los intervalos de confianza de la media de los tratamientos.
exp(confint.default(enfermedad_glm_1))## 2.5 % 97.5 %
## (Intercept) 0.003414506 58.115482
## hemo 1.249111045 2.461275
## bili 0.199751711 1.041384#Predecir
#Podemos predecir la probabilidad de encontrar un individuo en la N, con características como hemoglobina en 7 y bilirubina=3
pi.hat<- predict.glm(enfermedad_glm_1, data.frame(hemo=7, bili=3), type="response", se.fit=TRUE)
pi.hat## $fit
## 1
## 0.6828827
##
## $se.fit
## 1
## 0.1870155
##
## $residual.scale
## [1] 1pi.hat$fit## 1
## 0.6828827#La probabilidad de encontrar un individuo en la N es 68.29%Práctica
#En este ejemplo vamos a predecir si la presencia o ausencia de accidentes de automóvil (variable "presencia") en diversos tramos de carretera depende de las características de tramo (variable "orden", cuatro órdenes) y de la precipitación de lluvia (variable 'Precipitacion') utilizando regresión logística.
accidente <- read.table(url("http://www.uv.es/lejarza/actu/glm/trese.txt"), header = T, sep = "\t", dec = ",")
head(accidente)## Presencia montana Densidad usos Orden Precipitacion Sup.Tramo
## 1 1 -0.70 0.00 1.55 1 808 140.00
## 2 1 -0.70 0.00 2.28 2 1110 953.85
## 3 0 -0.70 0.05 0.18 1 754 41.17
## 4 0 -0.69 0.22 -4.65 1 421 109.42
## 5 1 0.08 0.19 -0.23 4 568 2009.94
## 6 0 0.08 0.00 -2.02 1 535 160.00str(accidente)## 'data.frame': 150 obs. of 7 variables:
## $ Presencia : int 1 1 0 0 1 0 0 0 0 1 ...
## $ montana : num -0.7 -0.7 -0.7 -0.69 0.08 0.08 0.08 0.08 0.08 0.08 ...
## $ Densidad : num 0 0 0.05 0.22 0.19 0 0.17 0.11 0 0.21 ...
## $ usos : num 1.55 2.28 0.18 -4.65 -0.23 -2.02 -0.52 0.42 0.21 0.41 ...
## $ Orden : num 1 2 1 1 4 1 2 1 1 3 ...
## $ Precipitacion: num 808 1110 754 421 568 535 570 690 668 556 ...
## $ Sup.Tramo : num 140 953.9 41.2 109.4 2009.9 ...plot(density(accidente$Precipitacion[accidente$Presencia == 0], adjust = 3), main = "", xlab = "Precipitacion", ylim = c(0, 0.004),lwd = 2)
lines(density(accidente$Precipitacion[accidente$Presencia == 1], adjust = 3), main = "", xlab = "Precipitacion", ylim = c(0, 0.004),lwd = 2, lty = 3)
#en estos gráficos así creados lo punteado se refiere a valores 1(si) de presencia de accidentes. En precipitación casi no hay diferencias , en tramossplot(density(accidente$Orden[accidente$Presencia == 0], adjust = 3), main = "", xlab = "Orden", lwd = 2)
lines(density(accidente$Orden[accidente$Presencia == 1], adjust = 3),main = "", xlab = "", ylab = "", lty = 3, lwd = 2) 
#Observando los datos .No parece haber una diferencia muy clara en los valores de precipitación entre tramos con y sin presencia de accidentes. Sin embargo, la presencia de accidentes parece estar asociada, a primera vista, a tramos con número (orden) mayor #Ajustando el modelo
accidente$Orden <- as.factor(accidente$Orden)
glm1 <- glm(Presencia ~ Precipitacion + Orden, data = accidente, family = binomial)
anova(glm1, test = "Chi") ## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Presencia
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 149 207.94
## Precipitacion 1 6.288 148 201.66 0.01215 *
## Orden 3 44.646 145 157.01 1.1e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Para observar los valores estimados de los parámetros así como las devianzas y el AIC, es mas detallado con summary.
summary(glm1)##
## Call:
## glm(formula = Presencia ~ Precipitacion + Orden, family = binomial,
## data = accidente)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0596 -0.8620 -0.2077 0.8227 1.9165
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.214e+00 1.663e+00 -3.135 0.001717 **
## Precipitacion 6.674e-03 2.448e-03 2.727 0.006393 **
## Orden2 1.870e+00 4.795e-01 3.900 9.61e-05 ***
## Orden3 3.950e+00 1.067e+00 3.700 0.000215 ***
## Orden4 1.699e+01 1.455e+03 0.012 0.990687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 207.94 on 149 degrees of freedom
## Residual deviance: 157.01 on 145 degrees of freedom
## AIC: 167.01
##
## Number of Fisher Scoring iterations: 14#interpretación.
#Ahora bien, no todos los coeficientes del factor tramo ("Orden"") son significativos. En principio, el Intercept, que resume el efecto del nivel de "Orden1"" sobre la presencia de accidentes, es significativo y negativo. Esto indicaría que en niveles de "Orden1"" la probabilidad de presencia de accidentes es menor que en el resto de niveles.
# Los niveles "Orden2"" y "Orden3"" también son significativos pero positivos, lo que indicaría que en estos tramos/niveles se incrementa presencia de accidentes.
#Por último, el nivel de "Orden4"" no es significativo, lo que indicaría que el valor positivo del coeficiente no es significativamente distinto de cero y, por tanto, tiene un efecto nulo sobre la variable respuesta. #Por último, es interesante saber qué proporción de la varianza explica el modelo.
(207.94-157.01)/207.94*100## [1] 24.49264# modelo explica el 24,49 % de la variabilidad.#podemo utilizar la función plot
par(mfrow=c(2,2))
plot(glm1)## Warning: not plotting observations with leverage one:
## 5
## Warning: not plotting observations with leverage one:
## 5
#Para ver residuos influyentes, pero no para la comprobación de los modelos#glm utilizando familia poisson
data("polyps", package = "HSAUR")
head(polyps)## number treat age
## 1 63 placebo 20
## 2 2 drug 16
## 3 28 placebo 18
## 4 17 drug 22
## 5 61 placebo 13
## 6 1 drug 23polyps_glm_1 <- glm(number ~ treat + age, data = polyps, family = poisson())
summary(polyps_glm_1)##
## Call:
## glm(formula = number ~ treat + age, family = poisson(), data = polyps)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.2212 -3.0536 -0.1802 1.4459 5.8301
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.529024 0.146872 30.84 < 2e-16 ***
## treatdrug -1.359083 0.117643 -11.55 < 2e-16 ***
## age -0.038830 0.005955 -6.52 7.02e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 378.66 on 19 degrees of freedom
## Residual deviance: 179.54 on 17 degrees of freedom
## AIC: 273.88
##
## Number of Fisher Scoring iterations: 5polyps_glm_2 <- glm(number ~ treat + age, data = polyps, family = quasipoisson())
summary(polyps_glm_2)##
## Call:
## glm(formula = number ~ treat + age, family = quasipoisson(),
## data = polyps)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.2212 -3.0536 -0.1802 1.4459 5.8301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.52902 0.48106 9.415 3.72e-08 ***
## treatdrug -1.35908 0.38533 -3.527 0.00259 **
## age -0.03883 0.01951 -1.991 0.06284 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 10.72805)
##
## Null deviance: 378.66 on 19 degrees of freedom
## Residual deviance: 179.54 on 17 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5plot(polyps_glm_2)
