MANOVA con variables de diferentes conjuntos de datos

Question

He realizado varios experimentos en los que tengo una variable independiente subyacente (especie de árbol, IV). Cada uno de estos experimentos me dio una variable dependiente (VD), como el pH de la corteza, la rugosidad o la capacidad de retención de agua. Ahora quiero realizar un MANOVA para ver si las especies de árboles difieren en las distintas variables dependientes. Mi análisis se realiza en R .

Por lo tanto, mi modelo es así:

pH + rugosity + water-holding capacity + [...] ~ tree species

donde tengo por especie de árbol ...

• 3 mediciones del pH de la corteza.
• 9 mediciones de la rugosidad de la corteza.
• 4 mediciones del espesor de la corteza,
• 5 mediciones de la capacidad de retención de agua,
• 5 mediciones de la retención de agua.

Sin embargo, a diferencia de la mayoría de los ejemplos que he encontrado sobre cómo hacer un MANOVA (es decir aquí , aquí , aquí ), mis datos provienen de diferentes mediciones y de diferentes individuos. Ahora, sólo he encontrado este hilo en el que se habla de tamaños de muestra desiguales, pero esto se refiere sólo a los tamaños de muestra dentro del factor explicativo.

Mi pregunta:

Todas mis variables dependientes tienen diferentes tamaños de muestra. ¿Sería apropiado un MANOVA para este tipo de datos? ¿Puedo simplemente ignorar los diferentes tamaños de las variables? ¿Existe una forma alternativa de hacerlo o más bien una prueba estadística alternativa? ¿Importa el pequeño tamaño de mi muestra?

EDIT: Lo que realmente quiero averiguar

En realidad sólo quiero realizar una prueba estadística que me diga si tengo un patrón subyacente. Entonces, ¿las especies de árboles son diferentes en cuanto a las variables dependientes? Al final quiero ser capaz de decir, si algunas especies tienen un cierto conjunto de rasgos diferentes de otras especies.

Ejemplo de datos:

Mis datos son así:

> manova_df
    # A tibble: 45 x 6
   tree_species rugosity bark_mm    pH   whc   ret
   <fct>           <dbl>   <int> <dbl> <dbl> <dbl>
 1 AS              2.36        8  6.49  295. 119. 
 2 AS              1.45        8  6.83  222. 105. 
 3 AS              3.13        9  5.8   291. 181. 
 4 AS              2.38        8 NA     314. 214. 
 5 AS              4.39        7 NA     613. 317. 
 6 AS              2.21       NA NA      NA   NA  
 7 AS              0.810      NA NA      NA   NA  
 8 AS              1.58       NA NA      NA   NA  
 9 AS              0.934      NA NA      NA   NA  
10 BU              3.34        6  7.22  189.  74.9
# ... with 35 more rows

El NA s provienen del hecho de que tengo diferentes tamaños de muestra, pero tenía que meter todas las variables en una sola data.frame . Así que simplemente uní todas las columnas de las diferentes observaciones. Esto significa que las observaciones individuales de las diferentes VD no están unidas y el orden dentro de los niveles de las especies del árbol es totalmente arbitrario.

Mi análisis es bastante sencillo:

mano_mod = manova(cbind(pH, bark_mm, rugosity, whc, ret) ~ tree_species, data = manova_df)
> summary(mano_mod)
             Df Pillai approx F num Df den Df    Pr(>F)    
tree_species  4 2.4207    3.372     20     44 0.0003836 ***
Residuals    12                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1  

Conjunto de datos completo:

structure(list(tree_species = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L), .Label = c("AS", "BU", "CL", 
"MB", "PR"), class = "factor"), rugosity = c(2.36, 1.45, 3.13, 
2.38, 4.39, 2.21, 0.81, 1.58, 0.93, 3.34, 5.06, 0, 0.77, 12.64, 
4.1, 0.8, 1.03, 0.84, 6.49, 9.09, 5.96, 5.32, 8.41, 15.29, 9.91, 
7.65, 2.13, 9.43, 10.14, 13.24, 10.26, 9.81, 12.34, 17.23, 16.63, 
8.82, 1.68, 0.7, 0.82, 2.43, 0, 0.76, 0.77, 0, 1), bark_mm = c(8L, 
8L, 9L, 8L, 7L, NA, NA, NA, NA, 6L, 8L, 8L, 7L, 9L, NA, NA, NA, 
NA, 9L, 9L, 8L, 10L, 9L, NA, NA, NA, NA, 5L, 9L, 9L, 8L, 4L, 
NA, NA, NA, NA, 5L, 5L, 5L, 6L, NA, NA, NA, NA, NA), pH = c(6.49, 
6.83, 5.8, NA, NA, NA, NA, NA, NA, 7.22, 7.11, 7.72, 7.29, NA, 
NA, NA, NA, NA, 7.39, 7.18, 7.3, 7.3, NA, NA, NA, NA, NA, 6.76, 
6.55, 6.24, NA, NA, NA, NA, NA, NA, 5.76, 6.59, 5.44, NA, NA, 
NA, NA, NA, NA), whc = c(295.2, 222.4, 290.6, 314.3, 613.4, NA, 
NA, NA, NA, 189.4, 248.2, 336.8, 330.1, 427.8, NA, NA, NA, NA, 
236, 492.6, 549.3, 330.1, 370.7, NA, NA, NA, NA, 430, 142.2, 
372.4, 260, 176.1, 680, 215, NA, NA, 333.8, 320.6, 282.4, 322.9, 
576.7, NA, NA, NA, NA), ret = c(118.9, 104.9, 180.6, 214.5, 317.3, 
NA, NA, NA, NA, 74.9, 95.7, 127.3, 150.1, 327.3, NA, NA, NA, 
NA, 80.8, 176.7, 255.7, 142.6, 236.6, NA, NA, NA, NA, 148.4, 
32.4, 244.2, 66.8, 76.4, 246.1, 73.6, NA, NA, 111.2, 151.3, 102.1, 
200.6, 258.1, NA, NA, NA, NA)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -45L))

(Si algo no está claro, pregúntelo).

Answer 1

Una opción es pensar en esto como un problema de datos perdidos. Tienes cantidades variables de información para cada caso (diferentes especies de árboles) que has clasificado como una de las cinco especies diferentes. Si puede cumplir con la suposición de que los datos faltan al azar, podría imputar.

Sin embargo, antes de seguir adelante, debo decirle que voy a argumentar en contra de un MANOVA para estos datos. Se trata de un conjunto de datos pequeño para que falten tantas puntuaciones como ya puedo ver en las primeras 10 filas. De hecho, si mira sus grados de libertad, sólo tiene un total de 17 casos con datos completos. El MANOVA va a suprimir la lista, por lo que se pierden todos los casos no completos. Ahora mismo su prueba inferencial se basa en menos de la mitad de sus datos y eso me parece subóptimo.

Mi opinión es que usted preferiría hacer inferencias utilizando la mayor cantidad de información posible, lo que va a ser limitado si está utilizando un MANOVA sin algún tipo de imputación. Incluso si se hace la imputación, Finch (2016) encontró que las tasas de falta superiores al 40%, incluso bajo la condición de que fueran faltas al azar (es decir, MAR), pueden causar una inflación indeseable del error de tipo I. Su estudio de simulación se basó en dos grupos con 50 observaciones por grupo (es decir, $N = 100$ ). Es probable que su modelo pueda tolerar incluso menos datos perdidos, dado que la muestra total es más pequeña ( $N_{max} = 45)$ y un mayor número de grupos ( $k=5$ ).

Si no es posible incluir las puntuaciones observadas para los casos a los que les faltan valores en sus diferentes medidas, yo recomendaría considerar un enfoque en el que usted impute sus datos usando tanta información como tenga sobre sus casos (es decir, si hay variables adicionales a las que tenga acceso que sean completamente observadas en todos los casos), y luego realice una serie de ANOVAs de una vía. Compruebe el mice y Amelia paquetes en R para las opciones de imputación.

Otra razón para dejar de lado el modelo MANOVA en su totalidad, aparte de los datos que faltan, es que puede no ser un ajuste conceptual tan bueno. A menos que tenga una hipótesis sólida sobre las diferencias de grupo en la ubicación de las distribuciones mutlivariadas (es decir, sus centros), ganará poco realizando un MANOVA en primer lugar. Existe el mito de que la realización de un MANOVA le protegerá de algún modo contra la inflación de tipo I que resulta de las comparaciones múltiples, algo que realmente no lo hace todo bien .

Mi opinión es que lo más probable es que le interese saber en qué dimensiones se diferencian sus distintas especies. De todos modos, iba a tener que realizar ANOVAs unidireccionales de seguimiento para averiguarlo, así que ¿por qué no eliminar al intermediario del MANOVA, que no hace nada para controlar el error de tipo I de estos ANOVAs unidireccionales separados?

ACTUALIZACIÓN:

Creo que el gls de Heteroskedastic Jim es una solución inteligente si quiere mantener esto en un modelo. Pero requiere un modelo complejo para un problema simple, y sin que se especifiquen efectos aleatorios más complejos en el modelo, no estoy seguro de lo que el enfoque "compra" por encima del enfoque ANOVA más simple con comparaciones múltiples. Para mayor claridad (utilizando los datos que ha proporcionado), esto es lo que obtendría con un flujo de trabajo ANOVA para rugosity . Nota a efectos de comparación con el gls ejemplo anterior no estoy imputando ningún valor aquí.

fit.lm1<-lm(rugosity~tree_species, data=dat)
fit.lm1.aov<-aov(fit.lm1) 

En primer lugar, y lo que es más importante, la prueba general ómnibus es significativa:

             Df Sum Sq Mean Sq F value   Pr(>F)    
tree_species  4  763.1  190.77   23.49 4.78e-10 ***
Residuals    40  324.9    8.12                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Y ahora los coeficientes del modelo:

Call:
lm(formula = rugosity ~ tree_species, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.6756 -1.8456 -0.1467  0.9244  9.4644 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)       2.138      0.950   2.250 0.029992 *  
tree_speciesBU    1.038      1.343   0.772 0.444374    
tree_speciesCL    5.668      1.343   4.219 0.000137 ***
tree_speciesMB    9.851      1.343   7.333 6.48e-09 ***
tree_speciesPR   -1.231      1.343  -0.916 0.364958    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.85 on 40 degrees of freedom
Multiple R-squared:  0.7014,    Adjusted R-squared:  0.6715 
F-statistic: 23.49 on 4 and 40 DF,  p-value: 4.775e-10

Pero para dar un poco más de sentido a las diferencias entre sus 5 especies, un conjunto de $p$ Los valores basados en comparaciones múltiples podrían ayudar. Aquí estoy utilizando la corrección más conservadora de Bonferonni para reducir la inflación de errores de tipo I, pero existen alternativas.

pairwise.t.test(dat$rugosity, dat$tree_species, p.adj = "bonferroni")

Que vuelve:

Pairwise comparisons using t tests with pooled SD 

data:  dat$rugosity and dat$tree_species 

   AS      BU      CL      MB     
BU 1.0000  -       -       -      
CL 0.0014  0.0135  -       -      
MB 6.5e-08 7.7e-07 0.0341  -      
PR 1.0000  0.9903  7.7e-05 3.6e-09

P value adjustment method: bonferroni 

Esa es su variable más completa. Entonces, ¿qué sucede con una DV diferente? Bueno, con bark_mm se obtendría la siguiente salida:

fit.lm2<-lm(bark_mm~tree_species, data=dat)
fit.lm2.aov<-aov(fit.lm2)
summary(fit.lm2.aov)
             Df Sum Sq Mean Sq F value  Pr(>F)   
tree_species  4  34.01   8.502   5.056 0.00603 **
Residuals    19  31.95   1.682                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
21 observations deleted due to missingness

summary(fit.lm2)

Call:
lm(formula = bark_mm ~ tree_species, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0000 -0.3375  0.0000  0.8125  2.0000 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      8.0000     0.5799  13.795 2.38e-11 ***
tree_speciesBU  -0.4000     0.8201  -0.488  0.63133    
tree_speciesCL   1.0000     0.8201   1.219  0.23765    
tree_speciesMB  -1.0000     0.8201  -1.219  0.23765    
tree_speciesPR  -2.7500     0.8699  -3.161  0.00514 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.297 on 19 degrees of freedom
  (21 observations deleted due to missingness)
Multiple R-squared:  0.5156,    Adjusted R-squared:  0.4136 
F-statistic: 5.056 on 4 and 19 DF,  p-value: 0.006027

pairwise.t.test(dat$bark_mm, dat$tree_species, p.adj = "bonferroni")

    Pairwise comparisons using t tests with pooled SD 

data:  dat$bark_mm and dat$tree_species 

   AS     BU     CL     MB    
BU 1.0000 -      -      -     
CL 1.0000 1.0000 -      -     
MB 1.0000 1.0000 0.2473 -     
PR 0.0514 0.1414 0.0038 0.5865

P value adjustment method: bonferroni 

Ahora compare esto con el gls salida del modelo:

Coefficients:
                                  Value Std.Error   t-value p-value
(Intercept)                      8.0000   0.57993 13.794673  0.0000
tree_speciesBU                  -0.4000   0.82015 -0.487715  0.6267
tree_speciesCL                   1.0000   0.82015  1.219288  0.2252
tree_speciesMB                  -1.0000   0.82015 -1.219288  0.2252
tree_speciesPR                  -2.7500   0.86990 -3.161279  0.0020
measurepH                       -1.6267   0.61771 -2.633383  0.0096
measureret                     179.2400  37.30808  4.804321  0.0000
measurerugosity                 -5.8622   1.11299 -5.267085  0.0000
measurewhc                     339.1800  64.50607  5.258110  0.0000
tree_speciesBU:measurepH         1.3617   0.86708  1.570412  0.1191
tree_speciesCL:measurepH        -0.0808   0.86708 -0.093225  0.9259
tree_speciesMB:measurepH         1.1433   0.87357  1.308800  0.1932
tree_speciesPR:measurepH         2.3067   0.92044  2.506045  0.0136
tree_speciesBU:measureret      -31.7800  52.76159 -0.602332  0.5481
tree_speciesCL:measureret       -9.7600  52.76159 -0.184983  0.8536
tree_speciesMB:measureret      -59.3971  48.84873 -1.215940  0.2265
tree_speciesPR:measureret      -19.8300  52.76239 -0.375836  0.7077
tree_speciesBU:measurerugosity   1.4378   1.57401  0.913450  0.3629
tree_speciesCL:measurerugosity   4.6678   1.57401  2.965536  0.0037
tree_speciesMB:measurerugosity  10.8511   1.57401  6.893937  0.0000
tree_speciesPR:measurerugosity   1.5189   1.60049  0.949012  0.3446
tree_speciesBU:measurewhc      -40.3200  91.22536 -0.441982  0.6593
tree_speciesCL:measurewhc       47.5600  91.22536  0.521346  0.6031
tree_speciesMB:measurewhc      -21.0800  84.45884 -0.249589  0.8034
tree_speciesPR:measurewhc       22.8500  91.22582  0.250477  0.8027

¿Dónde está su bark_mm ¿variable en el modelo? Bueno, puede que no sea obvio inicialmente, pero bark_mm es en realidad su referencia measure . Así que tendrá que ajustar el valor de referencia en su gls modelo para obtener la comparación que desea. Esto no es una condena del enfoque, sólo una pequeña molestia, y una razón por la que creo que el flujo de trabajo de ANOVA es más sencillo para conseguir lo que se quiere al final. Y si el error de tipo I es una preocupación, los ANOVAs todavía pueden ser utilizados responsablemente sin inflar el error de tipo I, con las correcciones apropiadas.

Además, no olvide que tiene que sacar sus comparaciones por pares de las diferentes especies del gls Si se mantiene ese enfoque (que es ciertamente factible - y como un aparte el pairwise.t.test() está lejos de ser la única opción para las comparaciones de medias post hoc).

Referencia(s)

Finch, W. H. (2016). Datos perdidos e imputación múltiple en el contexto del análisis multivariante de la varianza. Revista de Educación Experimental , 84, 356-372. doi: 10.1080/00220973.2015.1011594

Answer 2

Creo que se puede tratar el problema como una regresión multinivel en lugar de un MANOVA. El primer paso será hacer que los datos sean aún más largos:

dat$ID <- 1:nrow(dat)
dat.l <- tidyr::gather(dat, measure, value, rugosity:ret)
dat.l <- na.omit(dat.l)
dat.l <- dat.l[order(dat.l$ID), ]
head(dat.l)
# # A tibble: 6 x 4
#   tree_species    ID measure   value
#   <fct>        <int> <chr>     <dbl>
# 1 AS               1 rugosity   2.36
# 2 AS               1 bark_mm    8   
# 3 AS               1 pH         6.49
# 4 AS               1 whc      295.  
# 5 AS               1 ret      119.  
# 6 AS               2 rugosity   1.45

Ahora cada fila pertenece a una medición específica, tenemos un identificador para el tipo de medición y la especie de árbol.

Un modelo de referencia razonable sería aquel en el que permitiéramos que el value variable para diferir por medida con respecto a la media y la varianza, ya que las medidas parecen estar en escalas muy diferentes:

library(nlme)

fit.0 <- gls(
  value ~ measure, dat.l, # ~ 1 | ID,
  weights = varIdent(form = ~ 1 | measure))

Un segundo paso sería incluir la especie arbórea como predictor:

fit.1 <- gls(
  value ~ measure + tree_species, dat.l, # ~ 1 | ID,
  weights = varIdent(form = ~ 1 | measure))

A continuación, podemos comparar ambos modelos mediante las técnicas habituales de comparación de modelos:

anova(fit.0, fit.1)
#       Model df      AIC      BIC    logLik   Test  L.Ratio p-value
# fit.0     1 10 1054.608 1083.660 -517.3037                        
# fit.1     2 14 1035.080 1075.333 -503.5402 1 vs 2 27.52702  <.0001

Hay una advertencia de la que puede librarse utilizando ML en lugar de REML, algo así como anova(update(fit.0, method = "ML"), update(fit.1, method = "ML")) . Pero si nos centramos en los resultados que tenemos, una prueba de ratio de verosimilitud sugiere que el modelo más complicado se ajusta mejor a los datos. Y el AIC sugiere que el modelo más complicado tiene un mejor rendimiento predictivo fuera de la muestra. Así que incluir el tipo de especie puede ayudarnos a entender/predecir mejor los datos más allá de entender la medida específica que estamos viendo.

Un modelo aún más complicado sería permitir la interacción entre especies y medidas:

fit.2 <- gls(
  value ~ tree_species * measure, dat.l, # ~ 1 | ID,
  weights = varIdent(form = ~ 1 | measure))

En este modelo, ajustamos 16 parámetros adicionales ya que hay 5 medidas y 5 especies. Podemos volver a utilizar el anova() para comparar los tres modelos. Este modelo tendrá tantos coeficientes que podría ser difícil de entender, pero si los métodos de comparación sugieren que este es un modelo mejor, parece que las especies difieren en las formas de las diferentes medidas. Es probable que esto ocurra ya que las medidas están en escalas muy diferentes, incluso si es el mismo árbol el que es consistentemente diferente a través de las medidas. Probablemente tendrá que graficar el modelo o dividirlo para entender mejor lo que está sucediendo.