# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
DATOS <- read.table("~/desktop/curso R 2012/Crecimiento.csv", header = T, sep = ",")
attach(DATOS)
names(DATOS)
## [1] "Mes" "Planta" "Altura"
source("~/desktop/curso r 2012/funciones.r")
# Una tabulación para visualizar los datos
tapply(Altura, list(Planta, Mes), mean)
## 1 2 3 4 5 6 7 8 9
## 1 24.0 32.0 39.0 42.5 48.0 54.5 61.0 65.0 72.0
## 2 22.5 30.5 40.5 45.0 51.0 58.5 64.0 72.0 78.0
## 3 22.5 28.0 36.5 41.0 47.5 55.0 61.0 68.0 76.0
## 4 24.0 31.5 39.5 44.5 51.0 56.0 59.5 64.0 67.0
## 5 24.5 31.5 37.0 42.5 48.0 54.0 58.0 63.0 65.5
## 6 23.0 30.0 35.5 41.0 48.0 51.5 56.5 63.5 69.5
## 7 22.5 28.5 36.0 43.5 47.0 53.5 59.5 67.5 73.5
## 8 23.5 30.5 38.0 41.0 48.5 55.0 59.5 66.5 73.0
## 9 20.0 27.5 33.0 39.0 43.5 49.0 54.5 59.5 65.0
## 10 25.5 32.5 39.5 47.0 53.0 58.5 63.0 69.5 76.0
## 11 24.5 31.0 40.5 46.0 51.5 57.0 62.5 69.5 76.0
## 12 24.0 29.0 39.0 44.0 50.5 57.0 61.5 68.0 73.5
## 13 23.5 30.5 36.5 42.0 47.0 55.0 59.0 65.5 73.0
## 14 21.5 30.5 37.0 42.5 48.0 52.5 58.5 63.0 69.5
## 15 25.0 32.0 38.5 44.0 51.0 59.0 66.0 75.5 86.0
## 16 21.5 28.5 34.0 39.5 45.0 51.0 58.0 64.5 72.5
## 17 31.0 38.0 48.0 54.0 60.0 62.0 66.5 75.5 84.0
## 18 27.5 32.5 36.0 43.0 49.5 52.5 56.0 61.0 64.0
## 19 30.0 37.0 45.0 51.0 58.0 63.0 67.5 74.5 81.0
## 20 26.0 32.0 40.5 45.5 52.5 55.5 62.5 69.5 74.0
## 21 26.0 32.5 39.5 44.0 48.0 54.5 58.0 66.0 73.0
## 22 28.5 35.5 41.5 47.5 54.0 59.5 63.5 71.0 78.5
## 23 26.5 34.5 42.0 48.5 55.5 62.0 68.0 76.5 85.0
## 24 27.5 33.5 41.0 45.0 50.5 56.0 62.5 71.0 78.0
## 25 22.5 27.0 33.5 38.5 41.0 49.0 56.0 64.0 68.0
## 26 22.0 26.5 32.5 38.5 43.5 50.5 56.5 63.5 68.5
## 27 23.5 29.0 35.5 40.0 45.0 50.0 56.5 63.0 67.5
## 28 22.5 29.5 36.5 42.0 45.0 55.0 61.0 68.0 72.0
## 29 27.5 34.5 42.0 47.5 53.0 63.0 72.0 79.0 85.5
## 30 23.5 28.0 33.0 37.0 38.5 48.0 52.5 62.0 64.5
## 31 24.5 30.0 38.5 42.0 47.5 54.0 62.5 71.5 77.0
## 32 24.5 31.5 40.5 46.5 51.5 61.5 68.5 77.5 84.5
## 33 24.5 32.0 39.0 45.0 51.0 55.5 61.5 69.0 75.5
## 34 24.0 32.5 40.0 48.0 54.5 61.5 68.0 74.5 81.0
## 35 24.0 31.5 38.5 44.0 51.5 57.5 64.0 72.5 79.0
## 36 24.5 32.5 39.5 44.5 52.5 56.5 62.0 67.5 72.5
## 37 24.5 32.0 38.5 44.0 50.0 56.0 63.5 69.5 76.0
## 38 25.5 33.0 41.5 47.0 55.5 60.5 66.5 77.0 82.0
## 39 25.5 32.0 39.0 45.5 51.0 57.5 63.5 72.0 78.5
## 40 25.0 31.0 36.5 43.0 50.5 55.0 62.5 69.0 75.5
## 41 26.5 30.5 33.0 39.0 43.5 49.5 56.5 61.0 65.0
## 42 24.0 32.0 39.0 44.5 50.0 56.0 63.0 67.5 74.0
## 43 24.5 31.0 37.5 43.5 48.0 56.0 62.5 66.5 70.5
## 44 27.0 34.5 42.0 48.5 53.0 60.0 67.0 73.0 76.0
## 45 31.0 39.0 47.5 51.0 57.0 64.0 71.0 77.0 80.5
## 46 27.0 33.5 40.0 46.5 53.0 60.0 66.5 72.5 80.0
## 47 29.5 37.0 46.0 52.5 60.0 67.5 76.0 81.5 88.0
## 48 28.5 36.0 42.5 49.0 55.0 63.5 72.0 78.5 85.5
# las dimensiones de la base datos
dim(DATOS)
## [1] 432 3
# 432 obervaciones en 48 individuos, 9 observaciones por individuo
432/48
## [1] 9
# Graficamos medias vs varianzas, vemos que la varianza crece con la media
VAR <- tapply(Altura, Mes, var)
PROM <- tapply(Altura, Mes, mean)
par(mfrow = c(1, 1))
plot(VAR ~ PROM)
# Con GLM notamos una falta de homocedasticidad
M <- glm(Altura ~ Mes)
par(mfrow = c(2, 2))
plot(M)
# Cambiamos a error Gamma para corregir la falta de homocedasticidad, pero
# ahora hay sesgo en los residuales, se sobreestima en los extremos, y
# subestima en la parte central
M <- glm(Altura ~ Mes, Gamma)
par(mfrow = c(2, 2))
plot(M)
# Cambiamos de función de liga, de inverse a log, pero siguen habiendo
# sesgo en las estimaciones
M <- glm(Altura ~ Mes, Gamma("log"))
par(mfrow = c(2, 2))
plot(M)
# cambiamos a otra función de liga, sqrt, parece mejorar, pero no mucho,
# el sesgo persiste
M <- glm(Altura ~ Mes, Gamma("sqrt"))
par(mfrow = c(2, 2))
plot(M)
# Podemos incorporar el hecho de que hay 48 individuos, y mejora otro poco
# y podría ser aceptable, pero en la tabla de ANOVA ...
M <- glm(Altura ~ Planta + Mes, Gamma("sqrt"))
par(mfrow = c(2, 2))
plot(M)
# ... cabe preguntarse si realmente tenemos tantas replicas para la
# variable mes, 429!!!! solo hay 9 meses de registros
anova(M, test = "Chi")
## Analysis of Deviance Table
##
## Model: Gamma, link: sqrt
##
## Response: Altura
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 431 51.8
## Planta 1 0.4 430 51.4 3.5e-12 ***
## Mes 1 48.2 429 3.2 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Se esta pseudoreplicando el experimento ya que solo hay registros de 48
# plantas, cada una con nueve mediciones, estas mediciones no son
# independientes dentro de cada individuos, son medidas repetidas y esta
# medida repetida debe de ser considerado un factor aleatorio para evitar
# inflar los grados de libertad del modelo. En este ejemplo lo colocaremos
# como un factor aleatorio, un evento de medición, en el que se midio cada
# planta.
library(nlme)
MEM <- lme(Altura ~ Mes, random = ~1 | Mes/Planta)
summary(MEM)
## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 2517 2537 -1253
##
## Random effects:
## Formula: ~1 | Mes
## (Intercept)
## StdDev: 0.0002183
##
## Formula: ~1 | Planta %in% Mes
## (Intercept) Residual
## StdDev: 4.392 0.05602
##
## Fixed effects: Altura ~ Mes
## Value Std.Error DF t-value p-value
## (Intercept) 19.36 0.4605 423 42.03 0
## Mes 6.21 0.0818 7 75.88 0
## Correlation:
## (Intr)
## Mes -0.889
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -0.0345704 -0.0073606 -0.0005668 0.0075352 0.0382584
##
## Number of Observations: 432
## Number of Groups:
## Mes Planta %in% Mes
## 9 432
# Los residuales no parecen seguir el supuesto de la homocedasticidad,
# pero por un momento nos brincamos este paso para explorar la tabla de
# ANOVA de este modelo
plot(MEM)
# Ahora si hay solo siete grados de libertad,
anova(MEM)
## numDF denDF F-value p-value
## (Intercept) 1 423 56898 <.0001
## Mes 1 7 5757 <.0001
# y para cumplir los supuestos de los residuales dentro de cada mes,
# usamos una función de potencia para ponderar la varianza en cada mes
MEM <- lme(Altura ~ Mes, random = ~1 | Mes/Planta, weights = varPower())
# lo cual genera que los residuales se distribuyen de acuerdo a los
# supuestos de los modelos
plot(MEM)
anova(MEM)
## numDF denDF F-value p-value
## (Intercept) 1 423 46125 <.0001
## Mes 1 7 5981 <.0001
# y estos son las betas del predictor lineal
summary(MEM)$coefficients$fixed
## (Intercept) Mes
## 19.21 6.24
# Ahora comapre colocando mes del otro lado de la 'caña'
MEM2 <- lme(Altura ~ Mes, random = ~Mes | Planta)
# Aquí los residuales ya se distribuyen de forma adecuada si ponderar la
# varianza entre meses, aquí mes también es continuo en el factor
# aleatorio, en el ejemplo anterior fue categórico simplemente por
# colocarlo al lado derecho de la 'caña'
plot(MEM2)
# Ahora si hay solo siete grados de libertad,
anova(MEM2)
## numDF denDF F-value p-value
## (Intercept) 1 383 3292 <.0001
## Mes 1 383 4552 <.0001
dim(DATOS)
## [1] 432 3
432 - 48 - 1 #48 individuos
## [1] 383
# Estas son dos aproximaciones diferentes en el componente aleatorio del
# modelo, pero si comparos el resultado de las betas del factor fijo,
# vemos que estas no cambian
summary(MEM)$coefficients$fixed
## (Intercept) Mes
## 19.21 6.24
summary(MEM2)$coefficients$fixed
## (Intercept) Mes
## 19.36 6.21
# Con la primera aproximación podemos explorar si la identidad de las
# plantas tiene efecto en el modelo. Removemos del << Planta >> de la
# parte aleatoria del modelo,
MEM3 <- lme(Altura ~ Mes, random = ~1 | Mes, weights = varPower())
# Y la prueba de ANOVA nos indica que no hay efecto significativo de este
# factor aleatorio, su efecto en la varianza no es distinguible de cero
anova(MEM3, MEM)
## Model df AIC BIC logLik Test L.Ratio p-value
## MEM3 1 5 2455 2475 -1222
## MEM 2 6 2457 2481 -1222 1 vs 2 0.1524 0.6962
anova(MEM3)
## numDF denDF F-value p-value
## (Intercept) 1 423 52254 <.0001
## Mes 1 7 6993 <.0001
summary(MEM3)
## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 2455 2475 -1222
##
## Random effects:
## Formula: ~1 | Mes
## (Intercept) Residual
## StdDev: 0.007522 0.1769
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 0.8125
## Fixed effects: Altura ~ Mes
## Value Std.Error DF t-value p-value
## (Intercept) 19.178 0.3211 423 59.72 0
## Mes 6.248 0.0747 7 83.63 0
## Correlation:
## (Intr)
## Mes -0.826
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.78650 -0.58129 -0.04422 0.60807 2.96985
##
## Number of Observations: 432
## Number of Groups: 9
# Para graficar
SEC <- seq(1, 9, length = 30)
summary(MEM3)
## Linear mixed-effects model fit by REML
## Data: NULL
## AIC BIC logLik
## 2455 2475 -1222
##
## Random effects:
## Formula: ~1 | Mes
## (Intercept) Residual
## StdDev: 0.007522 0.1769
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~fitted(.)
## Parameter estimates:
## power
## 0.8125
## Fixed effects: Altura ~ Mes
## Value Std.Error DF t-value p-value
## (Intercept) 19.178 0.3211 423 59.72 0
## Mes 6.248 0.0747 7 83.63 0
## Correlation:
## (Intr)
## Mes -0.826
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.78650 -0.58129 -0.04422 0.60807 2.96985
##
## Number of Observations: 432
## Number of Groups: 9
par(mfrow = c(1, 1), mar = c(4, 4, 1, 1))
plot(Altura ~ Mes, type = "p", ylab = "Altura", xlab = "Mes", cex = 0.8, pch = 19)
# y usamos los valores de las betas del predictor lineal
summary(MEM3)$coefficients$fixed
## (Intercept) Mes
## 19.178 6.248
PRED <- 19.17764 + 6.248044 * SEC
lines(SEC, PRED, col = "red", lwd = 3)
# Esta linea de la predicción no se corresponde exactamente con la de los
# promedios de cada mes!, pero si se parece mucho
MEDIAS <- tapply(Altura, Mes, mean)
lines(1:9, MEDIAS, col = "blue", lwd = 3)
### Otro ejemplo con la base de datos precargada en el biblioteca << nlme
### >>
library(nlme)
data(Orthodont)
attach(Orthodont)
# Exploramos los datos
tapply(distance, list(Subject, age, Sex), mean)
## , , Male
##
## 8 10 12 14
## M16 22.0 21.5 23.5 25.0
## M05 20.0 23.5 22.5 26.0
## M02 21.5 22.5 23.0 26.5
## M11 23.0 23.0 23.5 25.0
## M07 22.0 22.0 24.5 26.5
## M08 24.0 21.5 24.5 25.5
## M03 23.0 22.5 24.0 27.5
## M12 21.5 23.5 24.0 28.0
## M13 17.0 24.5 26.0 29.5
## M14 22.5 25.5 25.5 26.0
## M09 23.0 20.5 31.0 26.0
## M15 23.0 24.5 26.0 30.0
## M06 24.5 25.5 27.0 28.5
## M04 25.5 27.5 26.5 27.0
## M01 26.0 25.0 29.0 31.0
## M10 27.5 28.0 31.0 31.5
## F10 NA NA NA NA
## F09 NA NA NA NA
## F06 NA NA NA NA
## F01 NA NA NA NA
## F05 NA NA NA NA
## F07 NA NA NA NA
## F02 NA NA NA NA
## F08 NA NA NA NA
## F03 NA NA NA NA
## F04 NA NA NA NA
## F11 NA NA NA NA
##
## , , Female
##
## 8 10 12 14
## M16 NA NA NA NA
## M05 NA NA NA NA
## M02 NA NA NA NA
## M11 NA NA NA NA
## M07 NA NA NA NA
## M08 NA NA NA NA
## M03 NA NA NA NA
## M12 NA NA NA NA
## M13 NA NA NA NA
## M14 NA NA NA NA
## M09 NA NA NA NA
## M15 NA NA NA NA
## M06 NA NA NA NA
## M04 NA NA NA NA
## M01 NA NA NA NA
## M10 NA NA NA NA
## F10 16.5 19.0 19.0 19.5
## F09 20.0 21.0 22.0 21.5
## F06 20.0 21.0 21.0 22.5
## F01 21.0 20.0 21.5 23.0
## F05 21.5 23.0 22.5 23.5
## F07 21.5 22.5 23.0 25.0
## F02 21.0 21.5 24.0 25.5
## F08 23.0 23.0 23.5 24.0
## F03 20.5 24.0 24.5 26.0
## F04 23.5 24.5 25.0 26.5
## F11 24.5 25.0 28.0 28.0
# Hay 16 niños y 11 niñas medidos en las edades de 8, 9 10, 11, 12 13, y
# 14 años
# Con un modelo lineal obtenemos una buena distribución de los residuales
M <- glm(distance ~ age * Sex)
par(mfrow = c(2, 2))
plot(M)
# Pero otra vez cabe la pregunta, ¿realmente tenemos tantos grados de
# libertad? La respuesta es no!
anova(M, test = "Chi")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: distance
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 107 918
## age 1 235.4 106 682 1.1e-11 ***
## Sex 1 140.5 105 542 1.5e-07 ***
## age:Sex 1 12.1 104 530 0.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Usamos el modelo mixto
MEM <- lme(distance ~ age * Sex, random = ~age | Subject)
# los residuales se ven bien
plot(MEM)
# Note el cambio en los grados de libertad
anova(MEM)
## numDF denDF F-value p-value
## (Intercept) 1 79 4036 <.0001
## age 1 79 99 <.0001
## Sex 1 25 8 0.0090
## age:Sex 1 79 5 0.0264
# Como hay mas de un termino en la parte fija, es conveniente usar el
# argumento << type = 'marginal' >> para que los cálculos se hagan con la
# suma de cuadrados marginales, similar a la simplificación de los modelo,
# y no con el modelo final. Compare el cambio en los valores de F y en la
# significancia
anova(MEM)
## numDF denDF F-value p-value
## (Intercept) 1 79 4036 <.0001
## age 1 79 99 <.0001
## Sex 1 25 8 0.0090
## age:Sex 1 79 5 0.0264
anova(MEM, type = "marginal")
## numDF denDF F-value p-value
## (Intercept) 1 79 257.39 <.0001
## age 1 79 83.19 <.0001
## Sex 1 25 0.42 0.5237
## age:Sex 1 79 5.12 0.0264
# Usamos la función predict de forma similar que lo usamos en los GLM y LM
EDAD <- seq(8, 14, length = 30)
# Y tenemos una predicción para cada sujeto de cada sexo!
ND <- data.frame(age = EDAD, Sex = "Female", Subject = "F11")
predict(MEM, ND)
## F11 F11 F11 F11 F11 F11 F11 F11 F11 F11 F11 F11
## 24.22 24.34 24.46 24.58 24.70 24.82 24.94 25.06 25.18 25.30 25.42 25.54
## F11 F11 F11 F11 F11 F11 F11 F11 F11 F11 F11 F11
## 25.66 25.78 25.90 26.02 26.14 26.26 26.38 26.50 26.62 26.74 26.86 26.98
## F11 F11 F11 F11 F11 F11
## 27.10 27.22 27.34 27.46 27.58 27.70
## attr(,"label")
## [1] "Predicted values"
ND <- data.frame(age = EDAD, Sex = "Female", Subject = "F01")
predict(MEM, ND)
## F01 F01 F01 F01 F01 F01 F01 F01 F01 F01 F01 F01
## 20.21 20.30 20.39 20.48 20.57 20.66 20.75 20.84 20.93 21.02 21.11 21.20
## F01 F01 F01 F01 F01 F01 F01 F01 F01 F01 F01 F01
## 21.29 21.38 21.47 21.56 21.65 21.74 21.83 21.92 22.01 22.10 22.19 22.28
## F01 F01 F01 F01 F01 F01
## 22.37 22.46 22.55 22.64 22.73 22.82
## attr(,"label")
## [1] "Predicted values"
# Entonces podemos graficar una linea para cada sujeto observado
par(mfrow = c(1, 1))
plot(age, distance, ylab = "Distancia entre premolares", xlab = "Edad (años)",
type = "n")
# Para las niñas
ESTOS <- levels(Subject)
for (i in ESTOS[17:27]) {
ND <- data.frame(age = EDAD, Sex = "Female", Subject = i)
lines(EDAD, predict(MEM, ND), col = "orchid")
}
# y para los niños
for (i in ESTOS[1:16]) {
ND <- data.frame(age = EDAD, Sex = "Male", Subject = i)
lines(EDAD, predict(MEM, ND), col = "navy")
}
# Sin embargo puede resulta mas útil graficar los datos, puntos, y la
# linea promedio pro sexo
plot(age, distance, ylab = "Distancia entre premolares", xlab = "Edad (años)",
type = "n")
ESTOS <- levels(Subject)
PRED <- numeric(30)
# Para las niñas
for (i in ESTOS[17:27]) {
EDAD <- seq(0, 14, length = 30)
ND <- data.frame(age = EDAD, Sex = "Female", Subject = i)
PRED <- PRED + as.vector(predict(MEM, ND))
}
lines(EDAD, PRED/11, , col = "orchid", lwd = 3)
points(age[Sex == "Female"], distance[Sex == "Female"], pch = 19, col = "orchid")
# Para los niños
PRED <- numeric(30)
for (i in ESTOS[1:16]) {
ND <- data.frame(age = EDAD, Sex = "Male", Subject = i)
PRED <- PRED + as.vector(predict(MEM, ND))
}
lines(EDAD, PRED/16, , col = "navy", lwd = 3)
points(age[Sex == "Male"] + 0.1, distance[Sex == "Male"], pch = 19, col = "navy")
# Note que los puntos correspondientes a los niños está ligeramente
# desfasados a la derecha simplemente para visualizarlos
# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Alternativamente podemos usar las betas del predictor lineal para pintar
# las mismas lineas
plot(age, distance, ylab = "Distancia entre premolares", xlab = "Edad (años)",
type = "n")
points(age[Sex == "Female"], distance[Sex == "Female"], pch = 19, col = "orchid")
points(age[Sex == "Male"] + 0.1, distance[Sex == "Male"], pch = 19, col = "navy")
# Note las betas que se multiplican por cero
summary(MEM)$coefficients$fixed
## (Intercept) age SexFemale age:SexFemale
## 16.3406 0.7844 1.0321 -0.3048
NIÑA <- 16.340625 + 0.784375 * EDAD + 1.032102 + -0.30483 * EDAD
NIÑO <- 16.340625 + 0.784375 * EDAD + 1.032102 * 0 + -0.30483 * EDAD * 0
lines(EDAD, NIÑA, col = "orchid", lwd = 3)
lines(EDAD, NIÑO, col = "navy", lwd = 3)
