Como modelar con heterogeneidad de varianza

En este ejemplo mostraré como usar los GLS para utilizarlos en vez de un ANOVA para ver diferencias en las medias de poblaciones con heterogeneidad de varianza. En resumen lo que se hace es incorporar la varianza en la ecuación.

Generación de base de datos de ejemplo:

Lo primero que haremos es generar una base de datos para mostrar como ejemplo:

library(tidyverse)
set.seed(1232)
pop = data.frame(group=c("A","B","C"),
                 mean=c(1,2,5),
                 sd=c(1,3,4))
d = do.call(rbind, rep(list(pop),13)) 
d$x = rnorm(nrow(d), d$mean, d$sd)

d <- d %>% select(group, x) %>% arrange(group)

La base de datos d, tiene una variable de interes x y una variable que los agrupa en poblaciones group, si graficamos un boxplot de esto, veremos que las medias son distintas, pero la varianza también, donde A tiene menos varianza, por lo que es imposible hacer un ANOVA.

ggplot(d, aes(x = group, y = x)) + geom_boxplot(aes(fill = group)) + theme_classic()

usando gls

Para usar los modelos gls hay que instalar y cargar la libreria nlme, para agregar la varianza en el modelo usamos los pesos con la funcion varIdent.

library(nlme)
mod.gls = gls(x ~ group, data=d, weights=varIdent(form= ~ 1 | group))

Al ver los resumenes de este modelo vemos lo siguiente:

summary(mod.gls)
## Generalized least squares fit by REML
##   Model: x ~ group 
##   Data: d 
##        AIC      BIC    logLik
##   177.1685 186.6696 -82.58426
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | group 
##  Parameter estimates:
##        A        B        C 
## 1.000000 2.444532 3.913382 
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 0.509768 0.2816667 1.809829  0.0787
## groupB      1.736692 0.7439273 2.334492  0.0253
## groupC      5.422838 1.1376880 4.766542  0.0000
## 
##  Correlation: 
##        (Intr) groupB
## groupB -0.379       
## groupC -0.248  0.094
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2911080 -0.5312418 -0.1639290  0.8286991  1.9881590 
## 
## Residual standard error: 1.015564 
## Degrees of freedom: 39 total; 36 residual

En este resultado vemos que se detectan diferencias significativas de B con A y de C con A, dados los P de B y C, pero ¿Como vemos si hay diferencias entre B y C? Para esto tenemos que crear un par de funciones y utilizar el paquete multicomp.

Para crear las funciones solo tienes que correr el código a continuación:

model.matrix.gls <- function(object, ...)
    model.matrix(terms(object), data = getData(object), ...)
model.frame.gls <- function(object, ...)
  model.frame(formula(object), data = getData(object), ...)
terms.gls <- function(object, ...)
  terms(model.frame(object),...)

y posteriormente utilizar estas funciones con el paquete multicomp.

library(multcomp)
mod.gls.mc = glht(mod.gls, linfct = mcp(group = "Tukey"))
summary(mod.gls.mc)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: gls(model = x ~ group, data = d, weights = varIdent(form = ~1 | 
##     group))
## 
## Linear Hypotheses:
##            Estimate Std. Error z value Pr(>|z|)    
## B - A == 0   1.7367     0.7439   2.334   0.0481 *  
## C - A == 0   5.4228     1.1377   4.767   <0.001 ***
## C - B == 0   3.6861     1.2996   2.836   0.0119 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Ahí por fin vemos que si hay diferencias significativas entre todos los grupos.