# Paquetes necesarios #### library(nlme) # install.packages("nlme") para instalar el paquete library(car) # install.packages("car") para instalar el paquete library(nlraa) # install.packages("nlraa") para instalar el paquete #### Simulación de diferentes varianzas entre grupos o niveles de un factor ##################### # Problema: # Crecimiento de dos especies (Sp1 y Sp2) en dos sitios (Norte y Sur). # -> variable respuesta: crecimiento en altura (cm/año). # -> dos factores con dos niveles cada uno: 2 x 2 = 4 tratamientos o grupos. # Población: # - el crecimiento de la Sp1 en el norte es 15 cm/año, # - en el sur el crecimiento es 7 cm/año menor que en el norte, # - la Sp2 crece 1 cm/año más que la Sp1, # - en el norte la varianza es 25 (σ = 5 cm/año), # - la varianza en el sur es 25 veces menor que en el norte (es decir, σ = 1 cm/año). # Modelo: # ci ~ N(µi, σi)indep. # µi = β0 + β1 * Sp2i + β2 * Si # σi = σ^2 * (Ni + δ1 * Si)^2 # con β0 = 15 (crecimiento de la Sp1 en el norte) # β1 = -7 (efecto del sitio) # β2 = 1 (efecto de la especie) # β3 = 0 (no hay efecto interacción) # σ = 5 (desvío estándar en el norte) # δ1 = 1/5 (cociente entre el desvío estándar en el norte y el del sur) # donde c = crecimiento (cm/año) # Sp2 = especie 2 (variable "dummy") # N = sitio norte (variable "dummy") # S = sitio sur (variable "dummy") # i = unidad experimental # Muestra: # - 20 observaciones en cada tratamiento (n = 20 x 4 = 80) # Paso 1.1: Simulación de datos ------------- # Generamos el vector sitio sitio <- rep(c("Norte","Sur"), each=40) # los primeros 40 valores para el sitio norte y 40 en el sur # Generamos el vector especie especie <- rep(rep(c("Sp1","Sp2"), each=20), 2) # 20 valores en Sp1 y 20 en Sp2, 20 en Sp1 y 20 en Sp2 # de manera que ambos factores estén cruzados grupos <- data.frame(especie, sitio) grupos table(grupos[,c(1,2)]) # grupos # Simulamos 80 valores de crecimiento de acuerdo al órden de los grupos (Sp1-N, Sp2-N, Sp1-S, Sp2-S) set.seed(23456) # para reproducir la aleatorización y1 <- rnorm(20, 15, 5) # Sp1 en el norte: 20 observaciones -> y1 ~ N(15;5) y2 <- rnorm(20, 16, 5) # Sp2 en el norte: 20 observaciones -> y2 ~ N(16;5) y3 <- rnorm(20, 8, 1) # Sp1 en el sur: 20 observaciones -> y3 ~ N( 8;1) y4 <- rnorm(20, 9, 1) # Sp2 en el sur: 20 observaciones -> y4 ~ N( 9;1) crecimiento <- c(y1, y2, y3, y4) crecimiento <- round(crecimiento, 2) # Asignamos los 4 grupos a los 80 crecimientos datos.sim <- data.frame(crecimiento, grupos) str(datos.sim) # 80 datos de 3 variables head(datos.sim) # Paso 1.2: Visualización de datos ------------ par(mfrow=c(1,1), cex=0.9) boxplot(crecimiento ~ especie + sitio, data=datos.sim, xlab="Grupo", ylab="Crecimiento (cm/año)", sep=" - ", col="grey", cex.lab=1.25) # heterocedasticidad entre sitios # Paso 1.3: Selección de la función de varianzas y ajuste del modelo ----------------- # Modelo lineal clásico (varianza cte) mlc <- gls(crecimiento ~ sitio * especie, # "sitio" y "especie" son los predictores del modelo (con "*" incluímos la interacción) data=datos.sim) # Modelo varIdent mvi <- gls(crecimiento ~ sitio * especie, weights = varIdent(form = ~ 1 | sitio), # varianza en función del sitio data=datos.sim) # Podemos chequear las matrices de varianzas y covarianzas de ambos modelos: vc.mvi <-round(var_cov(mvi),2) # matriz var-cov mvi vc.mvi[c(1:10), c(1:10)] # las primeras 10 observaciones corresponden al norte # y su varianza estimada es 27,3 (s=5,2) vc.mvi[c(71:80), c(71:80)] # las ultimas 10 observaciones al sur # y su varianza estimada es cercana a 1 (s≈1) vc.mlc <-round(var_cov(mlc),2) # matriz var-cov mlc vc.mlc[c(1:10), c(1:10)] # en el norte el modelo estima una varianza de 14,2 (s=3,8) vc.mlc[c(71:80), c(71:80)] # la misma que en el sur # Como los modelos está anidados podemos realizar el test de cociente de verosimilitudes (Caja 2, Tabla C2-1) anova(mlc, mvi) # y resulta que modelar la varianza mejora significativamente el ajuste # Residuales (Caja 2, Figura C2-1) # "mlc" residuos <- resid(mlc) # residuos "ordinarios" predichos <- fitted(mlc) # predichos par(mfrow=c(2,2), cex =0.7) plot(predichos, residuos, ylab="Residuos ordinarios", xlab="Predichos", main="Residuos vs Predichos") abline(0,0, col="red") plot(datos.sim$sitio, residuos, ylab="Residuos ordinarios", xlab="Sitio", main="Residuos vs Sitio") hist(residuos, ylab="Freceuncia", xlab="Residuos ordinarios", main="Histograma de residuos") qqnorm(residuos, ylab="Cuantiles muestrales", xlab="Cuantiles teóricos", main="Gráfico cuantil-cuantil") qqline(residuos) ks.test(residuos,"pnorm", mean(residuos), sd(residuos)) # test de Kolmogorov # p=0,027 -> falta de normalidad # "mvi" residuos.p <- resid(mvi, type="p") # residuos "Pearson" predichos <- fitted(mvi) # predichos par(mfrow=c(2,2), cex =0.7) plot(predichos, residuos.p, ylab="Residuos Pearson", main="Residuos vs Predichos") # la heterocedasticidad mejora notablemente abline(0,0, col="red") plot(datos.sim$sitio, residuos.p, ylab="Residuos Pearson", xlab="Sitio", main="Residuos vs Sitio") hist(residuos.p, ylab="Freceuncia", xlab="Residuos Pearson", main="Histograma de residuos") qqnorm(residuos.p, ylab="Cuantiles muestrales", xlab="Cuantiles teóricos", main="Gráfico cuantil-cuantil") qqline(residuos.p) ks.test(residuos.p,"pnorm", mean(residuos.p), sd(residuos.p)) # test de Kolmogorov # p=0,94 -> se corrige el problema de normalidad # Paso 1.4: Inferencias ----------------------- anova(mvi) # - no hay interacción entre especie y sitio # - diferencias entre sitios # - diferencias entre especies summary(mvi) # la función "summary()" devuelve el modelo ajustado # desde la cual extraemos las cantidades de interés: round(summary(mvi)$tTable[1, ],2) # la Sp1 en el norte tiene un crecimiento # estimado de 14,7 cm/año (b0) round(summary(mvi)$tTable[2, ],2) # el crecimiento estimado en el sur es 6,5 # cm/año menor que en el norte (b2) round(summary(mvi)$tTable[3, ],2) # el crecimiento estimado de la Sp2 es 0,9 # cm/año mayor que la Sp1 (b1) round(summary(mvi)$sigma,2) # en el norte el desvío estándar estimado es de 5,2 # cm/año (√varianza base), round((summary(mvi)$sigma)^2,2) # es decir, una varianza de 27,3 # con el siguiente código extraemos el vector de parámetros de varianza vp.var <- coef(mvi$modelStruct$varStruct, unconstrained=FALSE, allCoef=TRUE) round(vp.var,3) # y observamos que en el sur el desvío estándar es casi cinco veces menor # que en el norte (δ1 estimado = 0,19). Para obtener el desvío estándar en # el sur, multiplicamos el desvío del norte (√varianza base) por δ1: s.sur <- summary(mvi)$sigma * vp.var[2] round(s.sur,3) # aproximadamente 1 cm/año # Comparamos la estimación de parámetros de ambos modelos: # ...componente determinística (Caja 2, Tabla C2-2) compareCoefs(mvi, mlc) # estimación puntual y error estándar round(intervals(mvi)$coef,2) # IC modelo varIdent round(intervals(mlc)$coef,2) # IC modelo lineal clásico # ...componentes de varianza round(intervals(mvi)$sigma,2) # IC modelo varIdent round(intervals(mvi)$varStruct,2) round(intervals(mlc)$sigma,2) # IC modelo lineal clásico # - las estimaciones puntuales coinciden, # - el modelo lineal clásico estíma un desvío estándar común de 3,7 cm/año, # - el modelo varIdent estíma un desvío estandar de 5,5 cm/año en el norte y de 0,99 (5,2*0,19) cm/año en el sur, # - esto causa que la precisión de las estimaciones de β0 y β1 difiera entre modelos, # - y también los valores de F y sus valores-p asociados, tablas.anova <- cbind(anova(mvi)[,c(2:3)], anova(mlc)[,c(2:3)]) colnames(tablas.anova) <- c("F mvi", "valor-p mvi", "F mlc", "valor-p mlc") round(tablas.anova,3) # (Caja 2, Tabla C2-3) # - consecuentemente, con el modelo lineal clásico cometeríamos un error de tipo II al # concluir equivocadamente que no hay efecto de la especie.