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(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,