Convergence of results of nlme lme() based repeated measures ANOVA between R and InVivoStat

Question

My aim is to run repeated measures ANOVA. I am using AR1 covariance structure to model within-subject correlations in dependent variable as AR1 covariance structure based model fits my data better than Compound symmetry based model and equally compared to unconstrained covariance structure based model. I am using R package nlme for my analyses with following syntax:

model <- formula(Sah_1 ~ 1 + Pre_CVS.sahharoosi.eelistus + Gender + Grupp + HC.LC + Time.f + Gender:Grupp +
                     Gender:HC.LC + Gender:Time.f + Grupp:HC.LC + Grupp:Time.f + HC.LC:Time.f +
                     Gender:Grupp:HC.LC +
                     Gender:Grupp:Time.f +
                     Gender:HC.LC:Time.f +
                     Grupp:HC.LC:Time.f +
                     Gender:Grupp:HC.LC:Time.f
)


lme_data <- lme(model, random = ~1|grupimajutuse.ID, correlation=corAR1(form=~as.numeric(Time.f)|grupimajutuse.ID,fixed = FALSE), 
                na.action = (na.omit),
                data = SAH,
                method= "REML")
anova(lme_data)

Time.f is repeated factor (my measures are separated by constant time intervals);
Gender; HC.LC; Grupp are between-subjects treatment factors
Pre_CVS.sahharoosi.eelistus is between subjects co-variate.

My problem is non-convergence of my results obtained in R with results I obtain fitting identically specified model in InVivoStat. InVivoStat is free GUI based statistical package that runs R as its computational engine and InVivoStat documentations points out that repeated measures ANOVA is performed by using nlme package’s lme() function. In my case main effects and lower level interactions differ in F-statistics and p-values between R and InVivoStat, whereas higher order interactions return identical results. Here are the results I obtain myself from R:

                            numDF denDF  F-value p-value
(Intercept)                     1   120 413392.7  <.0001
Pre_CVS.sahharoosi.eelistus     1    39      4.9  0.0331
Gender                          1    39      0.0  0.9391
Grupp                           1    39      6.6  0.0140
HC.LC                           1    39      0.9  0.3422
Time.f                          3   120     16.2  <.0001
Gender:Grupp                    1    39      1.3  0.2683
Gender:HC.LC                    1    39      3.8  0.0582
Gender:Time.f                   3   120      0.6  0.6150
Grupp:HC.LC                     1    39      0.0  0.9358
Grupp:Time.f                    3   120      2.7  0.0516
HC.LC:Time.f                    3   120      2.6  0.0560
Gender:Grupp:HC.LC              1    39      0.0  0.9484
Gender:Grupp:Time.f             3   120      0.2  0.9155
Gender:HC.LC:Time.f             3   120      0.4  0.7464
Grupp:HC.LC:Time.f              3   120      0.4  0.7281
Gender:Grupp:HC.LC:Time.f       3   120      0.0  0.9869

Here are the results from InVivoStat:

Here is the outline of steps I have tried myself without success:

Changing the anova type from sequential to marginal and/or changing the REML estimator to ML estimator does not work as I lose the convergence even in higher order interaction terms in respect to F-statistics and p-values. In detail, changing the anova(lme_data) to anova(update(lme_data,method="ML")) does not help. Neither do anova(update(lme_data,method="ML"), type = "marginal") or anova(update(lme_data,method="REML"), type = "marginal") work.

Neither does forcing the variability of responses at different time periods to be equal help:

lme_data <- lme(model, random = ~1|grupimajutuse.ID, correlation=corAR1(form=~as.numeric(Time.f)|grupimajutuse.ID,fixed = FALSE), 
                na.action = (na.omit),
                data = SAH,
                method= "REML",
                weights = varFixed(~as.numeric(Time.f))
                ) 

Does anybody have an idea or suggestions based on my outlined syntax how to achieve the convergence of results between R and InVivoStats in repeated measures analysis. I.e., I am trying to guess the model specification used in InVivoStat to carry out previously described analysis. Unfortunately the source code of InVivoStat is not available so I can’t directly compare my syntax with it.

Thanks in advance,

Answer 1

The Repeated Measures Overall Test results presented by InVivoStat are (to use SAS notation) Type III or marginal.

To reproduce the results you first need to change R's default contrasts:

options(contrasts=c(unordered="contr.sum", ordered="contr.poly"))

The libraries InVivoStat uses are: multcomp, nlme, contrast and lsmeans.

This is the syntax of the AR1 model (obviously the variable names are different):

Model<-lme(model, random=~1|subject, correlation=corAR1(value=0.999, 
    form=~as.numeric(Time)|subject, fixed =FALSE), data=data,
    na.action = (na.omit), method = "REML")

To get the overall tests table you need:

anova(Model, type="marginal")

Hope this helps!