Equation for mixed model / residual-error structure

[Sarah Maulburg]

Hi All,

I am struggeling with putting together the correct equation for my model (see below). The study design: Outcome Variable is weimus (measuring the fatigue in patients with multiple sclerosis), subjects/patients (IDs) were randomized into two groups (intervention-and controll group) , we measured the outcome at four different timepoints (baseline, after 3 weeks, 3 months, 6 months) in each subject. I am using the mixed command to analyse the repeated measures data (change over time, change between groups and interaction between group * time). The within-subject covariance structure was not compound symmetric, which was one reasons to analyse the data with the mixed model instead of a repeated measures ANOVA.

Syntax:
mixed weimus group##time || id:, var noconst residuals(unstr, t(time)) reml

I checked the Stata Manual in order to find the correct equation for the model, but did not find out how to account for the unstructured residual-error covariance matrix.
What I thought of so far is:

WEIMuSij = β₀ + β₁ group_ij + β₂ time_ij + β₃ group x time_ij + u_j + E_ij for i = t0, t1, t2, t3 and j = 1,...,64 IDs (patients/subjects)

1. Would this equation be correct? (How) can I include the fitted residual-error structure?
2. Is it correct to say that I am using a mixed model or would it be more correct to say that it is a random intercept model?
3. What assumptions must be fullfilled?

Thank you so much in advance.

Best regards,
Sarah

[Joseph Coveney]

1. Would this equation be correct? (How) can I include the fitted residual-error structure?

With the noconstant option in mixed's random effects equation, you got rid of the u_j, and so the equation that you wrote isn't quite representative of the model that you're fitting. Representation of the residuals' unstructured covariance matrix seems right.

2. Is it correct to say that I am using a mixed model or would it be more correct to say that it is a random intercept model?

I would say neither. You're essentially fitting a MANOVA model using mixed instead of manova. To see that, run the code below and look at the output. (Start looking at the "Begin here" comment; the top part just creates a fictitious dataset for illustration.)

3. What assumptions must be fullfilled?

The assumptions that you'd find in the textbooks when they talk about MANOVA, including multivariate normality of the residuals.

version 15.1

clear *

set seed `=strreverse("1436721")'
local N 35

tempname Corr
matrix input `Corr' =  (1.000 0.750 0.500 0.250 \ ///
                        0.750 1.000 0.750 0.500 \ ///
                        0.500 0.750 1.000 0.750 \ ///
                        0.250 0.500 0.750 1.000)

// Control treatment group
quietly drawnorm weimus0 weimus3 weimus12 weimus24, ///
    means(50 50 50 50) sd(15 15 15 15) double corr(`Corr') n(`N')
generate byte group = 0
generate int pid = _n

tempfile tmpfil0
quietly save `tmpfil0'

// Experimental treatment group
quietly drawnorm weimus0 weimus3 weimus12 weimus24, ///
    means(50 40 40 40) sd(15 15 15 15) double corr(`Corr') n(`N') clear
generate byte group = 1
generate int pid = _n + `N'

append using `tmpfil0'

quietly reshape long weimus, i(pid) j(time)
summarize weimus, meanonly
quietly replace weimus = round(68 * (weimus - r(min)) / (r(max) - r(min)))

*
* Begin here
*
// -mixed-
mixed weimus i.group##i.time || pid: , reml dfmethod(kroger) ///
    noconstant residuals(unstructured, t(time)) ///
    nolrtest nolog

* Main effects of treatment
contrast r.group, small pveffects

* Main effects of time
contrast r.time, small pveffects

* Treatment × time interaction
contrast r.group#r.time, small pveffects

// -manova-
quietly reshape wide weimus, i(pid) j(time)

manova weimus0-weimus24 = group

matrix input Test = (-1 1 0)
matrix input Notest = (0.5 0.5 1)

matrix input Contrast = (1 -1 0 0 \ 1 0 -1 0 \ 1 0 0 -1)
matrix define Nocontrast = J(1, 4, 1)

* Main effects of treatment
manovatest , test(Test) ytransform(Nocontrast)

* Main effects of time
manovatest , test(Notest) ytransform(Contrast)

* Treatment × time interaction
manovatest , test(Test) ytransform(Contrast)

exit