Importancia del componente en un modelo lineal multivariante ("dentro del diseño")

Question

Estoy intentando probar la importancia del efecto "componente" en un modelo de regresión multivariante. No estoy seguro de cuál es la forma correcta. Usando R, he probado una forma con lm() y de otra manera con gls() y no dan resultados compatibles.

Tenga en cuenta que esta no es una pregunta sobre qué metodología es la correcta para analizar mis datos. Por cierto, se trata de datos simulados. Mi pregunta es sobre la comprensión en términos matemáticos de los procedimientos de R que utilizo.

El conjunto de datos:

> str(dat)
'data.frame':   31 obs. of  5 variables:
 $ group: Factor w/ 5 levels "1","5","2","3",..: 1 1 1 1 1 1 1 1 3 3 ...
 $ id   : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 3 ...
 $ x    : num  2.5 3 3 4 1.2 3.8 3.9 4 2.5 2.9 ...
 $ y    : num  2.6 3.8 3.9 3.8 1.6 5.2 1.3 3.6 4 3.2 ...
 $ z    : num  3.1 3.6 4.9 3.8 2.1 6 2.1 2.9 4.2 2.9 ...
> head(dat,10)
   group id   x   y   z
1      1  1 2.5 2.6 3.1
2      1  2 3.0 3.8 3.6
3      1  3 3.0 3.9 4.9
4      1  4 4.0 3.8 3.8
5      1  5 1.2 1.6 2.1
6      1  6 3.8 5.2 6.0
7      1  7 3.9 1.3 2.1
8      1  8 4.0 3.6 2.9
9      2  1 2.5 4.0 4.2
10     2  3 2.9 3.2 2.9

Convierto este conjunto de datos en "formato largo" para los gráficos (y posteriormente para gls() ):

dat$subject <- dat$group : dat$id
dat.long <- reshape(dat, direction="long", varying=list(3:5), 
    idvar="subject", v.names="value", timevar="component", times=c("x","y","z"))
dat.long$component <- factor(dat.long$component)

xyplot(value ~ component | group, data=dat.long, 
    pch=16, 
    strip = strip.custom(strip.names=TRUE,var.name="group" ), layout=c(5,1))

Cada individuo de cada grupo tiene $3$ medidas repetidas $x$ , $y$ , $z$ (Debería unir los puntos en el gráfico para ver las medidas repetidas).

Quiero ajustar un modelo MANOVA utilizando el grupo como factor y $(x,y,z)$ es la respuesta multivariante: $$\begin{pmatrix} x_{ij} \\ y_{ij} \\ z_{ij} \end{pmatrix} \sim {\cal N}_3\left( \begin{matrix} \mu_{i1} \\ \mu_{i2} \\ \mu_{i3} \end{matrix}, \Sigma \right), \quad i=1,\ldots,5$$ (por supuesto, podríamos utilizar la parametrización por defecto de R $\mu_{ik}=\mu_{1k} + \alpha_{ik}$ considerando group1 como el "intercepto" para cada respuesta pero prefiero "mi" parametrización).

Este modelo se ajusta de la siguiente manera utilizando lm() :

###  multivariate least-squares fitting  ###
mfit <- lm( cbind(x,y,z)~group, data=dat )

Creo que el modelo también se puede equipar con gls() de la siguiente manera (pero con un procedimiento de ajuste diferente) :

### generalized least-squares fitting  ###
library(nlme)
gfit <- gls(value ~ group*component, data=dat.long, correlation=corSymm(form= ~ 1 | subject))

Recordemos que subject = group:id es el identificador de los individuos. La página web correlation=corSymm(form= ~ 1 | subject) significa que las respuestas $x$ , $y$ , $z$ para cada individuo están correlacionados. Aquí corSymm significa una estructura de covarianza general, "no restringida" (denominada "no estructurada" en el lenguaje de SAS).

Para comprobar que mfit y gfit son equivalentes, podemos comprobar, por ejemplo, que podemos deducir los parámetros estimados de mfit a partir de los parámetros estimados de gfit y viceversa (por lo que los parámetros "medios" tienen exactamente los mismos valores ajustados):

> coef(mfit)
                  x          y          z
(Intercept)  3.1750  3.2250000  3.5625000
group5      -0.9500 -0.4750000  0.1125000
group2      -1.0750 -0.5678571 -0.2339286
group3      -0.7875 -0.1000000  0.1875000
group4      -0.3750  0.4000000 -0.0125000
> coef(gfit)
      (Intercept)            group5            group2            group3 
        3.1750000        -0.9500000        -1.0750000        -0.7875000 
           group4        componenty        componentz group5:componenty 
       -0.3750000         0.0500000         0.3875000         0.4750000 
group2:componenty group3:componenty group4:componenty group5:componentz 
        0.5071429         0.6875000         0.7750000         1.0625000 
group2:componentz group3:componentz group4:componentz 
        0.8410714         0.9750000         0.3625000 

Ahora quiero probar el "efecto componente". Hablando con rigor, escribiendo el modelo como $$\begin{pmatrix} x_{ij} \\ y_{ij} \\ z_{ij} \end{pmatrix} \sim {\cal N}_3\left( \begin{matrix} \mu_{i1} \\ \mu_{i2} \\ \mu_{i3} \end{matrix}, \Sigma \right),$$ Quiero probar la hipótesis $\boxed{H_0\colon \{\mu_{1i}=\mu_{2i}=\mu_{3i} \quad \forall i=1,2,3,4,5 \}}$ .

A continuación están mis intentos, un intento con gfit y dos intentos con mfit() :

###########################################
## testing significance of the component ##
###########################################

> ### with gfit  ###
> anova(gfit)
Denom. DF: 78 
                numDF  F-value p-value
(Intercept)         1 504.5226  <.0001
group               4   0.7797  0.5418
component           2  11.7073  <.0001
group:component     8   0.7978  0.6063
> 
> ### with mfit ###
> library(car)
> 
> # first attempt : 
> idata <- data.frame(component=c("x","y","z"))
> ( av.ok <- Anova(mfit, idata=idata, idesign=~component, type="III") )

Type III Repeated Measures MANOVA Tests: Pillai test statistic
                Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)      1   0.84396  140.625      1     26 5.449e-12 ***
group            4   0.10369    0.752      4     26    0.5658    
component        1   0.04913    0.646      2     25    0.5328    
group:component  4   0.22360    0.818      8     52    0.5901    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> # second attempt :
> linearHypothesis(mfit, "(Intercept) = 0", idata=idata, idesign=~component, iterms="component")

 Response transformation matrix:
  component1 component2
x          1          0
y          0          1
z         -1         -1

Sum of squares and products for the hypothesis:
           component1 component2
component1    1.20125    1.04625
component2    1.04625    0.91125

Sum of squares and products for error:
           component1 component2
component1   31.46179   14.67696
component2   14.67696   21.42304

Multivariate Tests: 
                 Df test stat  approx F num Df den Df  Pr(>F)
Pillai            1 0.0491253 0.6457903      2     25 0.53277
Wilks             1 0.9508747 0.6457903      2     25 0.53277
Hotelling-Lawley  1 0.0516632 0.6457903      2     25 0.53277
Roy               1 0.0516632 0.6457903      2     25 0.53277

Con anova(gfit) el componente es significativo, pero no con mis dos intentos usando mfit y el car paquete.

Sé que gls() utilizar un método de adaptación diferente al de lm() pero seguramente esta no es la causa de la diferencia.

Así que mis preguntas son :

• ¿He hecho algo mal?
• cuyo método prueba mi $H_0$ ¿hipótesis?
• qué es el $H_0$ ¿hipótesis de los otros métodos?

Y tengo una pregunta auxiliar: ¿cómo conseguir $\hat\Sigma$ con mfit y gfit ?

Actualización 1

A continuación se presenta un ejemplo reproducible que simula el conjunto de datos. Ahora creo que entiendo : ambos métodos ANOVA son correctos (el primero con anova(gfit) y el segundo con Anova(mfit, ...) , y arrojan resultados muy cercanos cuando se utiliza la suma de cuadrados de tipo II en Anova(mfit, ...) . Para el ejemplo anterior:

> Anova(mfit, idata=idata, idesign=~component, type="II")

Type II Repeated Measures MANOVA Tests: Pillai test statistic
                Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)      1   0.94691   463.72      1     26 < 2.2e-16 ***
group            4   0.10369     0.75      4     26 0.5657691    
component        1   0.50344    12.67      2     25 0.0001584 ***
group:component  4   0.22360     0.82      8     52 0.5900848    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

está muy cerca de

> anova(gfit)
Denom. DF: 78 
                numDF  F-value p-value
(Intercept)         1 504.5226  <.0001
group               4   0.7797  0.5418
component           2  11.7073  <.0001
group:component     8   0.7978  0.6063

A continuación se muestra el código reproducible con el muestreador de datos (simulo medidas repetidas no correlacionadas pero basta con incluir una matriz de covarianza en el rmvnorm() para simular medidas repetidas correlacionadas) :

library(mvtnorm)
library(nlme)
library(car)

# set data parameters 
I <- 5 # number of groups
J <- 16 # number of individuals per group
dat <- data.frame(
    group = gl(I,J),
    id = gl(J,1,I*J),
    x=NA, 
    y=NA, 
    z=NA
)
Mu <- c(1:I) # group means of components (assuming E(x)=E(y)=E(z) in each group)

# simulates data: 
for(i in 1:I){
    which.group.i <- which(dat$group==i)
    dat[which.group.i,c("x","y","z")] <- round(rmvnorm(n=J, mean=rep(Mu[i],3)), 1)
}

dat$subject <- droplevels( dat$group : dat$id )
dat.long <- reshape(dat, direction="long", varying=list(3:5), 
    idvar="subject", v.names="value", timevar="component", times=c("x","y","z"))
dat.long$component <- factor(dat.long$component)

# multivariate least-squares fitting 
mfit <- lm( cbind(x,y,z)~group, data=dat )

# gls fitting
dat.long$order.xyz <- as.numeric(dat.long$component)
gfit <- gls(value ~ group*component , data=dat.long, correlation=corSymm(form=  ~ order.xyz | subject)) 

# compares ANOVA : 
anova(gfit)
idata <- data.frame(component=c("x","y","z"))
Anova(mfit, idata=idata, idesign=~component, type="II")
Anova(mfit, idata=idata, idesign=~component, type="III")

Así que ahora me pregunto qué tipo de suma de cuadrados es el más apropiado para mi estudio real... pero esta es otra cuestión

Actualización 2

Sobre mi pregunta "cómo conseguir $\hat\Sigma$ " aquí está la respuesta para gls() :

> getVarCov(gfit)
Marginal variance covariance matrix
        [,1]    [,2]    [,3]
[1,] 0.92909 0.47739 0.24628
[2,] 0.47739 0.92909 0.53369
[3,] 0.24628 0.53369 0.92909
  Standard Deviations: 0.96389 0.96389 0.96389 

Eso demuestra que mfit y gfit no eran modelos equivalentes : gfit asume la misma varianza para los tres componentes.

Para ajustar una matriz de covarianza totalmente sin restricciones para las medidas repetidas, tenemos que escribir:

gfit2 <- gls(value ~ group*component , data=dat.long, 
    correlation=corSymm(form=  ~ 1 | subject), 
    weights=varIdent(form = ~1 | component))

> summary(gfit2)
Generalized least squares fit by REML
  Model: value ~ group * component 
  Data: dat.long 
       AIC      BIC    logLik
  264.0077 313.4986 -111.0038

Correlation Structure: General
 Formula: ~1 | subject 
 Parameter estimate(s):
 Correlation: 
  1     2    
2 0.529      
3 0.300 0.616
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | component 
 Parameter estimates:
       x        y        z 
1.000000 1.253534 1.169335 

....

Residual standard error: 0.8523997 

Pero aún no entiendo la matriz de covarianza extraída dada por getVarCov() (pero esto no es importante ya que obtenemos esta matriz con summary(gfit2) ):

   > getVarCov(gfit2)
    Error in t(S * sqrt(vars)) : 
      dims [product 9] do not match the length of object [0]
    > getVarCov(gfit2, individual="1:1")
    Marginal variance covariance matrix
            [,1]    [,2]    [,3]
    [1,] 0.72659 0.48164 0.25500
    [2,] 0.48164 1.14170 0.65562
    [3,] 0.25500 0.65562 0.99349
      Standard Deviations: 0.8524 1.0685 0.99674 
    > getVarCov(gfit2, individual="1:2")
    Marginal variance covariance matrix
            [,1]    [,2]    [,3]
    [1,] 1.14170 0.56319 0.27337
    [2,] 0.56319 0.99349 0.52302
    [3,] 0.27337 0.52302 0.72659
      Standard Deviations: 1.0685 0.99674 0.8524 

Lamentablemente, el anova(gfit2) mesa no está tan cerca de Anova(mfit, ..., type="II") como anova(gfit) :

> anova(gfit2)
Denom. DF: 78 
                numDF  F-value p-value
(Intercept)         1 498.1744  <.0001
group               4   1.0514  0.3864
component           2  13.1801  <.0001
group:component     8   0.8310  0.5780

> Anova(mfit, idata=idata, idesign=~component, type="II")

Type II Repeated Measures MANOVA Tests: Pillai test statistic
                Df test stat approx F num Df den Df    Pr(>F)    
(Intercept)      1   0.94691   463.72      1     26 < 2.2e-16 ***
group            4   0.10369     0.75      4     26 0.5657691    
component        1   0.50344    12.67      2     25 0.0001584 ***
group:component  4   0.22360     0.82      8     52 0.5900848    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

> anova(gfit)
Denom. DF: 78 
                numDF  F-value p-value
(Intercept)         1 504.5226  <.0001
group               4   0.7797  0.5418
component           2  11.7073  <.0001
group:component     8   0.7978  0.6063

Answer 1

Si entiendo bien, tienes datos de medidas repetidas: 3 réplicas de forma independiente en cada sujeto que ha sido aleatorizado a algún tratamiento.

Has estimado dos modelos de regresión, mfit que está llamando el modelo de regresión lineal multivariante con resultados x, y, z estimados conjuntamente en el factor de grupo. Entonces hay gfit que es el modelo de mínimos cuadrados generalizados con un intercepto aleatorio para medidas repetidas sobre el sujeto.

En mi experiencia con los análisis de medidas repetidas, nunca he visto que se utilice el modelo multivariante para probar este tipo de hipótesis. Me quedo con su convención de utilizar $\mu_i$ para representar el efecto del tratamiento (diferencia media) en el $i$ -a medida (sin embargo $\beta_i$ sería una opción más adecuada). Creo que la hipótesis que estás probando en el ANOVA de mfit es:

$\mu_i = \mu_j = \mu_k = 0$

La razón por la que esta inferencia es errónea es que impone la restricción de que el intercepto TIENE que ser 0 para el $i$ -a, $j$ -th, y $k$ -el grupo de tratamiento. Si el intercepto tiene un valor positivo, esto infla el error de tipo I (rechazando la hipótesis nula más a menudo cuando no debería). Además, este es un mal método para inferir los efectos del tratamiento porque los está tratando como resultados separados cuando le gustaría combinar sus resultados para un análisis más potente.

La opción GLS es el camino a seguir.

La especificación de su modelo es incorrecta y creo que debe eliminar el efecto fijo component porque su especificación como efecto aleatorio en la corstructura dará cuenta de las medidas repetidas dentro de los individuos. No debería estimar directamente las medias dentro de los componentes, porque la motivación de un modelo de efectos mixtos es dar cuenta de la variabilidad intrasujeto.

Sin embargo, la estimación de las diferencias medias dentro de cada período de medición (x, y, z) es un gran análisis exploratorio.