16/10, 2018
Repaso
Entren a la siguiente página
Entender la hipótesis estadística a generar basado en la hipótesis biológica
Pregunta 1
- ¿Cual de estas relaciones es la que uno esperaría ver al finalizar el experimento?
Pregunta 2
- ¿Cuál de estos códigos resultaría en el modelo que generaría este gráfico?
Pregunta 2 (cont)
library(lme4)
ChickPoissMM1 <- glmer(weight ~
Time + Diet + (1 | Chick),
family = poisson, data = ChickWeight)
ChickPoissMM2 <- glmer(weight ~
Time + Time:Diet + (1 | Chick),
family = poisson, data = ChickWeight)
ChickPoissMM3 <- glmer(weight ~
Time + Time * Diet + (1 | Chick),
family = poisson, data = ChickWeight)
Pregunta 3
- ¿Cuanto pesaria un pollo en el dia 10 en la dieta 2 basado en el modelo anterior?
| term | estimate | std.error | statistic | p.value | group |
|---|---|---|---|---|---|
| (Intercept) | 3.84 | 0.03 | 121.74 | 0 | fixed |
| Time | 0.07 | 0.00 | 63.72 | 0 | fixed |
| Time:Diet2 | 0.01 | 0.00 | 4.99 | 0 | fixed |
| Time:Diet3 | 0.02 | 0.00 | 12.59 | 0 | fixed |
| Time:Diet4 | 0.01 | 0.00 | 5.91 | 0 | fixed |
| sd_(Intercept).Chick | 0.21 | NA | NA | NA | Chick |
Pregunta 3 (cont)
\[y = beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2 + C_0 + C_1\]
Pregunta 3 (cont)
exp(3.83 + 10 * 0.07 + 0.02)
## [1] 94.63241
exp(3.83 + 10 * 0.07 + (10 * 0.02))
## [1] 113.2956
3.83 + 10 * 0.07 + 0.02
## [1] 4.55
3.83 + 10 * 0.07 + (10 * 0.02)
## [1] 4.73
Pregunta 3 discución
- Interpretación de los parámetros
| term | estimate | std.error | statistic | p.value | group |
|---|---|---|---|---|---|
| (Intercept) | 3.84 | 0.03 | 121.74 | 0 | fixed |
| Time | 0.07 | 0.00 | 63.72 | 0 | fixed |
| Time:Diet2 | 0.01 | 0.00 | 4.99 | 0 | fixed |
| Time:Diet3 | 0.02 | 0.00 | 12.59 | 0 | fixed |
| Time:Diet4 | 0.01 | 0.00 | 5.91 | 0 | fixed |
| sd_(Intercept).Chick | 0.21 | NA | NA | NA | Chick |
Pregunta 4
- ¿Cuál es el código que mejor representa este gráfico?
Pregunta 5
- ¿Cuantas pendientes e interceptos espero ver en este modelo?
Pregunta 6
- ¿Como interpreto ese modelo?
| term | estimate | std.error | statistic | p.value | group |
|---|---|---|---|---|---|
| (Intercept) | 3.67 | 0.05 | 79.47 | 0.00 | fixed |
| Time | 0.08 | 0.00 | 123.05 | 0.00 | fixed |
| Diet2 | 0.16 | 0.08 | 2.02 | 0.04 | fixed |
| Diet3 | 0.33 | 0.08 | 4.16 | 0.00 | fixed |
| Diet4 | 0.29 | 0.08 | 3.75 | 0.00 | fixed |
| sd_(Intercept).Chick | 0.20 | NA | NA | NA | Chick |
Criterios de información y crossvalidation
Pregunta 7
- Se realiza un experimento en un invernadero con 12 plantas, 6 de Quebec y 6 de Mississippi, de cada lugar de origen 3 se enfriaron y 3 no, luego se procedió a medir la captación de \(CO_2\) de cada planta a 95, 175, 250, 350, 500, 675 y 1000 ppm de \(CO_2\) ¿Que criterio de información debiera usar para elegir entre modelos que compiten?
Pregunta 7 (cont)
Pregunta 8
- Asumiendo que se cumplen los supuestos ¿Que modelo debo elegir entre estos dos?
ChickPoissMM <- glmer(weight ~
Time + Time:Diet + (1 | Chick),
family = poisson, data = ChickWeight)
ChickGammaMM <- glmer(weight ~
Time + Time:Diet + (1 | Chick),
family = Gamma, data = ChickWeight)
| logLik | AIC | BIC | deviance | df.residual | Family | link |
|---|---|---|---|---|---|---|
| -2638.24 | 5290.48 | 5321.00 | 24.39 | 571 | gamma | Inverse |
| -2645.09 | 5302.17 | 5328.33 | 1294.15 | 572 | Poisson | log |
Pregunta 9
- ¿Cuál de estos modelos tendra un AICc más bajo?
| (Intercept) | Agriculture | Catholic | Education | Examination | R^2 | df | logLik | |
|---|---|---|---|---|---|---|---|---|
| 8 | 86.225 | -0.203 | 0.145 | -1.072 | NA | 0.642 | 5 | -160.706 |
| 16 | 91.055 | -0.221 | 0.124 | -0.962 | -0.261 | 0.650 | 6 | -160.206 |
Control mas detallado de dredge en MuMin
Pregunta 10
¿Cuantas variables como máximo puedo tener en un modelo para explicar empleo en la base de datos longley?
| GNP.deflator | GNP | Unemployed | Armed.Forces | Population | Year | Employed | |
|---|---|---|---|---|---|---|---|
| 1947 | 83.0 | 234.29 | 235.6 | 159.0 | 107.61 | 1947 | 60.32 |
| 1948 | 88.5 | 259.43 | 232.5 | 145.6 | 108.63 | 1948 | 61.12 |
| 1949 | 88.2 | 258.05 | 368.2 | 161.6 | 109.77 | 1949 | 60.17 |
| 1950 | 89.5 | 284.60 | 335.1 | 165.0 | 110.93 | 1950 | 61.19 |
| 1951 | 96.2 | 328.98 | 209.9 | 309.9 | 112.08 | 1951 | 63.22 |
| 1952 | 98.1 | 347.00 | 193.2 | 359.4 | 113.27 | 1952 | 63.64 |
| 1953 | 99.0 | 365.38 | 187.0 | 354.7 | 115.09 | 1953 | 64.99 |
| 1954 | 100.0 | 363.11 | 357.8 | 335.0 | 116.22 | 1954 | 63.76 |
| 1955 | 101.2 | 397.47 | 290.4 | 304.8 | 117.39 | 1955 | 66.02 |
| 1956 | 104.6 | 419.18 | 282.2 | 285.7 | 118.73 | 1956 | 67.86 |
| 1957 | 108.4 | 442.77 | 293.6 | 279.8 | 120.44 | 1957 | 68.17 |
| 1958 | 110.8 | 444.55 | 468.1 | 263.7 | 121.95 | 1958 | 66.51 |
| 1959 | 112.6 | 482.70 | 381.3 | 255.2 | 123.37 | 1959 | 68.66 |
| 1960 | 114.2 | 502.60 | 393.1 | 251.4 | 125.37 | 1960 | 69.56 |
| 1961 | 115.7 | 518.17 | 480.6 | 257.2 | 127.85 | 1961 | 69.33 |
| 1962 | 116.9 | 554.89 | 400.7 | 282.7 | 130.08 | 1962 | 70.55 |
Como agrego eso a MuMIn
fit <- lm(Employed ~ ., data = longley)
options(na.action = "na.fail")
dd <- dredge(fit, m.lim = c(0,
floor(nrow(longley)/10)))
Como agrego eso a MuMIn
| (Intercept) | Armed.Forces | GNP | GNP.deflator | Population | Unemployed | Year | df | logLik | AICc | delta | weight |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 51.84 | NA | 0.03 | NA | NA | NA | NA | 3 | -14.90 | 37.81 | 0.00 | 0.98 |
| -1335.11 | NA | NA | NA | NA | NA | 0.72 | 3 | -19.30 | 46.60 | 8.79 | 0.01 |
| 33.19 | NA | NA | 0.32 | NA | NA | NA | 3 | -19.42 | 46.84 | 9.03 | 0.01 |
| 8.38 | NA | NA | NA | 0.48 | NA | NA | 3 | -21.84 | 51.68 | 13.87 | 0.00 |
| 59.29 | NA | NA | NA | NA | 0.02 | NA | 3 | -39.96 | 87.91 | 50.11 | 0.00 |
| 59.30 | 0.02 | NA | NA | NA | NA | NA | 3 | -40.41 | 88.82 | 51.01 | 0.00 |
| 65.32 | NA | NA | NA | NA | NA | NA | 2 | -42.29 | 89.49 | 51.69 | 0.00 |
Pregunta 11
- ¿Que problema tiene la siguiente base de datos para LMs o GLMs?
- Ejemplo base de datos
Cementde MuMIn
Como agrego eso a MuMIn
GlobalMod <- lm(y ~ X1 + X2 + X3 +
X4, data = Cement)
cor(Cement[, -1])
| X1 | X2 | X3 | X4 | |
|---|---|---|---|---|
| X1 | 1.0000000 | 0.2285795 | -0.8241338 | -0.2454451 |
| X2 | 0.2285795 | 1.0000000 | -0.1392424 | -0.9729550 |
| X3 | -0.8241338 | -0.1392424 | 1.0000000 | 0.0295370 |
| X4 | -0.2454451 | -0.9729550 | 0.0295370 | 1.0000000 |
Como agrego eso a MuMIn (cont)
nm <- colnames(Cement[-1])
smat <- abs(cor(Cement[, -1])) <=
0.7
smat[!lower.tri(smat)] <- NA
| X1 | X2 | X3 | X4 | |
|---|---|---|---|---|
| X1 | NA | NA | NA | NA |
| X2 | TRUE | NA | NA | NA |
| X3 | FALSE | TRUE | NA | NA |
| X4 | TRUE | FALSE | TRUE | NA |
Como agrego eso a MuMIn (cont)
options(na.action = "na.fail") Selected <- dredge(GlobalMod, subset = smat)
| (Intercept) | X1 | X2 | X3 | X4 | df | logLik | AICc | delta | weight |
|---|---|---|---|---|---|---|---|---|---|
| 52.58 | 1.47 | 0.66 | NA | NA | 4 | -28.16 | 69.31 | 0.00 | 0.84 |
| 103.10 | 1.44 | NA | NA | -0.61 | 4 | -29.82 | 72.63 | 3.32 | 0.16 |
| 131.28 | NA | NA | -1.20 | -0.72 | 4 | -35.37 | 83.74 | 14.43 | 0.00 |
| 72.07 | NA | 0.73 | -1.01 | NA | 4 | -40.96 | 94.93 | 25.62 | 0.00 |
| 117.57 | NA | NA | NA | -0.74 | 3 | -45.87 | 100.41 | 31.10 | 0.00 |
| 57.42 | NA | 0.79 | NA | NA | 3 | -46.04 | 100.74 | 31.42 | 0.00 |
| 81.48 | 1.87 | NA | NA | NA | 3 | -48.21 | 105.08 | 35.77 | 0.00 |
| 110.20 | NA | NA | -1.26 | NA | 3 | -50.98 | 110.63 | 41.31 | 0.00 |
| 95.42 | NA | NA | NA | NA | 2 | -53.17 | 111.54 | 42.22 | 0.00 |
Incluir además límite de numero de variables
options(na.action = "na.fail")
Selected <- dredge(GlobalMod, subset = smat,
m.lim = c(0, floor(nrow(Cement)/10)))
| (Intercept) | X1 | X2 | X3 | X4 | df | logLik | AICc | delta | weight |
|---|---|---|---|---|---|---|---|---|---|
| 117.57 | NA | NA | NA | -0.74 | 3 | -45.87 | 100.41 | 0.00 | 0.51 |
| 57.42 | NA | 0.79 | NA | NA | 3 | -46.04 | 100.74 | 0.33 | 0.43 |
| 81.48 | 1.87 | NA | NA | NA | 3 | -48.21 | 105.08 | 4.67 | 0.05 |
| 110.20 | NA | NA | -1.26 | NA | 3 | -50.98 | 110.63 | 10.22 | 0.00 |
| 95.42 | NA | NA | NA | NA | 2 | -53.17 | 111.54 | 11.13 | 0.00 |
Otra Forma de modelos Binomiales
Binomial
library(faraway)
data("orings")
logitmod <- glm(cbind(damage, 6 -
damage) ~ temp, family = binomial,
orings)
| null.deviance | df.null | logLik | AIC | BIC | deviance | df.residual |
|---|---|---|---|---|---|---|
| 38.9 | 22 | -14.84 | 33.67 | 35.95 | 16.91 | 21 |
| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 11.66 | 3.30 | 3.54 | 0 |
| temp | -0.22 | 0.05 | -4.07 | 0 |
pchisq(deviance(logitmod), df.residual(logitmod),
lower = FALSE)
## [1] 0.7164099
Predicción
Discusión de artículo
Mayor error tipo 1
- Error tipo 1
- Decir que hay un efecto cuando en realidad no lo hay
Elegir el modelo según propiedades de los datos y diagnósticos
- Usualmente hay heterocedasticidad
- a mayor valor de y, mayor variabilidad en residuales
- Puede haber sobredispersión
- Los datos pueden ser cero-inflados
Ejemplo de heterocedasticidad
Problemas de error tipo I
- Solo en bases de datos de n bajo (< 30)
- Solucionable con resampling