1
$\begingroup$

I'm trying to test the significance of the "component" effect in a multivariate regression model. I'm not sure what is the right way. Using R, I have tried a way with lm() and another way with gls(), and they don't yield compatible results.

Please note that this is not a question about which methodology is the right one to use to analyze my data. By the way these are simulated data. My question is about the understanding in mathematical terms of the R procedures I use.

The dataset:

> 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

I convert this dataset into "long format" for graphics (and later for 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))

enter image description here

Each individual of each group has $3$ repeated measures $x$,$y$,$z$ (I should join the points in the graphic to see the repeated measures).

I want to fit a MANOVA model using group as factor and $(x,y,z)$ is the multivariate response: $$ \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 $$ (of course we could use the default R parameterization $\mu_{ik}=\mu_{1k} + \alpha_{ik}$ by considering group1as the "intercept" for each response but I prefer "my" parameterization).

This model is fitted as follows using lm():

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

I think the model can also be fitted with gls() as follows (but with a different fitting procedure) :

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

Recall that subject = group:id is the identifier of the individuals. The correlation=corSymm(form= ~ 1 | subject) argument means that the responses $x$, $y$, $z$ for each individual are correlated. Here corSymm means a general, "unrestricted", covariance structure (termed as "unstructured" in SAS language).

To check that mfit and gfit are equivalent, we can check for instance that we can deduce the estimated parameters of mfit from the estimated parameters of gfitand vice-versa (so the "mean" parameters have exactly the same fitted values):

> 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 

Now I want to test the "component effect". Rigorously speaking, writing the model as $$ \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), $$ I want to test the hypothesis $\boxed{H_0\colon \{\mu_{1i}=\mu_{2i}=\mu_{3i} \quad \forall i=1,2,3,4,5 \}}$.

Below are my attempts, one attempt with gfit and two attempts with 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

With anova(gfit) the component is significant, but not with my two attempts using mfit and the car package.

I know that gls() use a different fitting method than lm() but this is surely not the cause of the difference.

So my questions are :

  • did I do something wrong ?
  • which method tests my $H_0$ hypothesis ?
  • what is the $H_0$ hypothesis of the other methods ?

And I have an auxiliary question: how to get $\hat\Sigma$ with mfit and gfit ?

Update 1

Below is a reproducible example which simulates the dataset. Now I think I understand : both ANOVA methods are correct (the first one with anova(gfit) and the second one with Anova(mfit, ...), and they yield very close results when using the type II sum of squares in Anova(mfit, ...). For the above example:

> 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 

is very close to

> 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

Below is the reproducible code with the data sampler (I simulate uncorrelated repeated measures but it suffices to include a covariance matrix in the rmvnorm() function to simulate correlated repeated measures) :

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")

So now I wonder which type of sum of squares is the more appropriate one for my real study... but this is another question

Update 2

About my question "how to get $\hat\Sigma$", here is the answer for 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 

That shows that mfitand gfit were not equivalent models: gfitassumes the same variance for the three components.

In order to fit a fully unrestricted covariance matrix for the repeated measures, we have to type:

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 

But yet I don't understand the extracted covariance matrix given by getVarCov() (but this is not important since we get this matrix with 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 

Unfortunately, the anova(gfit2) table is not as close to Anova(mfit, ..., type="II") as 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
$\endgroup$
-1
$\begingroup$

If I understand right, you have repeated measures data: 3 replications independently on each subject whose been randomized to some treatment.

You've estimated two regression models, mfit which you're calling the multivariate linear regression model with outcomes x, y, z jointly estimated on the group factor. Then there's gfit which is the generalized least squares model with a random intercept for repeated measures on subject.

In my experience with repeated measures analyses, I've never seen the multivariate model used to test this kind of hypothesis. I'll stick to your convention of using $\mu_i$ to represent the treatment effect (mean difference) on the subject's $i$-th measurement (however $\beta_i$ would be a more appropriate choice). I think the hypothesis you're actually testing in the ANOVA of mfit is:

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

The reason this inference is wrong is that it imposes the restriction that the intercept HAS to be 0 for the $i$-th, $j$-th, and $k$-th treatment group. If the intercept is positive valued, this inflates the type I error (rejecting the null hypothesis more often when you shouldn't). Furthermore, this is a bad method for inferring treatment effects because you're treating them as separate outcomes when you'd like to combine their results for a more powerful analysis.

The GLS option is the way to go.

Your model specification is wrong there and I think you need to remove the fixed effect component because your specification of it as a random effect in the corstructure will account for repeated measures within individuals. You shouldn't directly estimate the means within components, because the motivation for a mixed effect model is to account for intrasubject variability.

Estimating mean differences within each measurement period (x, y, z), though, is a great exploratory analysis.

$\endgroup$
  • $\begingroup$ Thanks for you effort. But your answer is a little confused. There's no random effect here, this is nothing but a classical multivariate regression with fixed effects (a "MANOVA"). The repeated measures are $(x,y,z)$ and the model assumes a multivariate normal distribution for this multivariate response, with a trivariate mean depending on the group. $\endgroup$ – Stéphane Laurent Jan 18 '13 at 21:57
  • $\begingroup$ You're correct the multivariate linear model has no random effect. The generalized least squares does. You specified a correlation structure. You may wish to read Longitudinal Data Analysis by Liang Heagerty Zeger Diggle. $\endgroup$ – AdamO Jan 19 '13 at 0:02
  • $\begingroup$ A correlation structure for the repeated measures is not a random effect. And you don't answer to any of my questions, except telling me I am wrong without mathematical explanations. $\endgroup$ – Stéphane Laurent Jan 19 '13 at 6:53
  • $\begingroup$ If your mathematical explanations are wrong, you need to point this out to facilitate the discussion. For instance, I am unsure whether you understand my point that your correlation structure is accounting for the influence of repeated measures and so adjusting for group is possibly inflating the type II error rate. About the similarities between correlation structures and random effects, they are both methods of handling repeated measures data and I believe they yield asymptotically equivalent inference on tests of fixed effects using MLE, Score tests, or Wald Tests. $\endgroup$ – AdamO Jan 19 '13 at 19:51
  • $\begingroup$ Using the subject as a random effect would yield a repeated measures model with a "compound symmetry" (exchangeable) covariance structure. Here the structure has no restriction. I have edited my post hoping it is clearer now. This discussion about the random effect is off-topic. Morevoer my question is a mathematical one and talking about "influence of repeated measures" and "inflation of type II error if I include group as factor" is off-topic too. The group is included as a factor in my question, I don't ask for any advice about this choice. $\endgroup$ – Stéphane Laurent Jan 19 '13 at 21:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.