What am I modeling with R vs. SPSS mixed model? [closed]

Question

I am trying to model a 3-way repeated measures experiment, FixedFactorA * FixedFactorB * Time[days]. There are no missing observations, but my groups (FactorA * FactorB) are unequal (close, but not completely balanced). From reading online, the best way to model a repeated measures experiment in which observation order matters (due to the response mean and variance changing in a time-dependent way) and for unequal groups is to use a mixed model and specify an appropriate covariance structure. However, I am new to the idea of mixed models and I am confused as to whether I am using the correct syntax to model what I am trying to model.

I would like to do a full factorial analysis, such that I could detect significant time * factor interactions. For example, for subjects with FactorA = 1, their responses over time might have a different slope and/or intercept than subjects with FactorA =2. I also want to be able to check whether certain combinations of FactorA and FactorB have significantly different responses over time (hence the full three-way interaction term).

From reading online, it seems like AR1 is a reasonable covariance structure for longitudinal-like data, so I decided to try that. Also, I saw that one is supposed to use ML if one plans to compare two different models, so I chose that approach in anticipation of needing to fine-tune the model. It is also my understanding that the goal is to minimize the AIC during model selection.

This is the code in the log for what I tried in SPSS (for long-form data), which yielded an AIC of 2471:

MIXED RESPONSE BY FactorA FactorB Day
  /CRITERIA=CIN(95) MXITER(100) MXSTEP(10) SCORING(1) SINGULAR(0.000000000001) HCONVERGE(0,
ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)
  /FIXED=FactorA FactorB Day FactorA*FactorB FactorA*Day FactorB*Day FactorA*FactorB*Day | SSTYPE(3)
  /METHOD=ML
  /PRINT=SOLUTION TESTCOV
  /REPEATED=Day | SUBJECT(Subject_ID) COVTYPE(AR1)

This is what I tried in R, which yielded an AIC of 2156:

    require(nlme)

    #output error fix: https://stats.stackexchange.com/questions/40647/lme-error-iteration-limit-reached

    ctrl <- lmeControl(opt='optim')

    fit1 <- lme(RESPONSE ~ Day*FactorA*FactorB, random = ~ Day|Subject_ID, control=ctrl,
        correlation=corAR1(form=~Day), data, method="ML")

    summary(fit1)

These are my questions:

1. The SPSS code above yielded a model with AIC = 2471, while the R code yielded a model with AIC = 2156. What is it about the codes that makes the models different?

2. From what I described above, are either of these models appropriate for what I am trying to test? If not, what would be a better way, and how would I do it in both programs to get the same results?

3. What is the appropriate way to compare models in both R and SPSS, such that the comparison tests are the same regardless of the software being used? For example, how do you do something like anova.lme() in SPSS?

Edits

Another thing to note is that I didn't dummy-code my factors. I don't know if this is a problem for either software, or if the built-in coding is different in SPSS vs R. I also don't know if this will be a problem for my three-way interaction term.

Also, when I say "factor", I mean an unchanging group or characteristic (like "sex").

Answer 1

Normally, CrossValidated is not for programming questions. That would be Stack Overflow, but it looks like your SPSS syntax does not include a random component nor have you specified a covariance structure, which would yield different AIC values.

For longitudinal analysis using multi-level modelling (also called hierarchical level modelling), you first should create an"unconditional model," such that

$Level 1:$

$y=\beta_0+\beta_1(time_{j})+e_{j}$

$Level 2:$

$\beta_0=\gamma_{00}+U_0$

$\beta_1=\gamma_{10}+U_1$

Where $Var(e_{j})=\sigma^2$, $Var(U_0)=\tau_{00})$, $Var(U_1)=\tau_{11}$, and $Cov(U_1,U_2)=\tau_{01}$, or in matrix notation:

\begin{equation*} \begin{bmatrix} \ \tau_{00} & \tau_{01} \\ \ \tau_{10} & \tau_{11} \end{bmatrix} \end{equation*}

Now, construct a "conditional model," with a control group ($group_{j2}$)

$Level 1:$

$y=\beta_0+\beta_1(time_{j1})+\beta_2(group_{j2})+e_{j}$

$Level 2:$

$\beta_0=\gamma_{00}+U_0$

$\beta_1=\gamma_{10}+U_1$

$\beta_2=\gamma_{20}+U_2$

Where $Var(e_{j})=\sigma^2$, but the $\tau$ matrix now has an added element, such that:

\begin{equation*} \begin{bmatrix} \ \tau_{00} & \tau_{01} &\tau_{02}\\ \ \tau_{10} & \tau_{11} & \tau_{12}\\ \ \tau_{20} & \tau_{21} & \tau_{22} \end{bmatrix} \end{equation*}

Further, you can add two time-variant predictors (which you call factors and I notate as $factor_ij$). Thus, observations over $time_j$ are nested in factors, such that:

$Level 1:$

$y=\beta_0+\beta_1(time_{j1})+\beta_2(group_{j2})+e_{ij}$

$Level 2:$

$\beta_0=\gamma_{00}+\gamma_{01}(factor_{ij})+U_0$

$\beta_1=\gamma_{10}+\gamma_{11}(factor_{ij})+U_1$

$\beta_2=\gamma_{20}+\gamma_{21}(factor_{ij})+U_2$

Where Where $Var(e_{j})=\sigma^2$, and the $\tau$ matrix is:

\begin{equation*} \begin{bmatrix} \ \tau_{00} & \tau_{01} &\tau_{02}\\ \ \tau_{10} & \tau_{11} & \tau_{12}\\ \ \tau_{20} & \tau_{21} & \tau_{22} \end{bmatrix} \end{equation*}

Then, repeat each of these models, but where $Var(e_{ij})=AR(1)$, with rows and columns equal to $j$

\begin{equation*} \begin{bmatrix} \ 1 & \rho & \dots & \rho^{j-1}\\ \ \rho & 1 & \dots & \rho^{j-2}\\ \ \vdots & \vdots &\ddots & \vdots\\ \ \rho^{j-1} & \dots & \dots & 1 \end{bmatrix} \end{equation*}

Each of these six models yields an AIC. You can compare using a $\chi^2$ difference test, where the AIC is your test statistic and the number of estimated parameters is your $df$.

For more information, see Raudenbush & Bryk (2002).