Partially Nested design and xtmixed

[Francesca Pesola]

Dear StataListers,

I am analysing data from an RCT with two study arms (control vs. therapy). Patients in the therapy arm received CBT therapy while those in the control group were on a waiting list (partially nested design). I am aware I should take into account the therapist effect as 2 therapists were involved in the study.

I am wondering whether -xtmixed- can handle this kind of design or whether there is an alternative command I could use in Stata 13.

Many thanks,
Francesca

[Scott Baldwin]

You can model partially nested designs with -mixed- in Stata. Basically, you fit a model with a random slope for treatment but suppress the overall random intercept. Something like:

mixed y tx || therapistid: tx, nocons

My colleagues and I discuss the model options and evaluate their performance in:

Baldwin, S. A., Bauer, D. J., Stice, E., & Rohde, P. (2011). Evaluating models for partially clustered designs. Psychological Methods, 16(2), 149–165. http://doi.org/10.1037/a0023464

Given that you only have two therapists, you don't have good options (estimating a variance with two therapists is tough) for a model. Hope this helps.

Best,
Scott

[Francesca Pesola]

Hi Scott,
Thank you for your reply and the reference. I have downloaded your paper.

In the case of partially nested design I can only have random slopes but no intercepts?

The other question is would you adjust for therapist effect in your primary analysis or would it be OK to do as part of a sensitivity analysis as therapist allocation was not taken into account when developing the design?

[Scott Baldwin]

You are fitting a random intercept for the treatment condition. Using the random slope method is needed to correctly specify the design matrix for the random effects in a program like -mixed-, which doesn't allow for the direct specification of a random intercept for some levels of a dichotomous variable but not others. The problem you have with just two therapists is that you likely have some non-independence in your data, but no good way to model it. Suppose the sensitivity analysis doesn't show a difference form the model without the therapists. Should we trust it? Likewise, suppose it shows a difference. Should we trust that? However, simply ignoring the therapists isn't a good idea either. Some recommend fitting fixed effects for therapists in your case. We show in our paper that the fixed effects approach has problems in some situations. Anyhow, no good answers from my perspective.

Best,
Scott

[Francesca Pesola]

Hi Scott,

Thanks again for your reply. After reading your paper, I feel that I should account for the clustering in the experimental group in my main analysis to avoid any bias; I set the mixed command as follows:

qui sum bY,d
g Ycentred = bY -(r(mean))
mixed fY arm Ycentred Centre covs1 cov2 cov3 || clust: arm, nocons vce(robust)

where clust codes the therapist in the experimental group (coded: 1 and 2) and it is set to missing (.) in the control group; arm codes the study arm; Ycentred are the baseline scores centred around the mean. I receive a message from Stata saying

note: arm omitted because of collinearity AND the coefficient for arm is not estimated but set to 0.

Should arm not be missing (.) in the control group?

A more general question is: is requesting robust SEs - when you have data collected at 2 time points - not a viable option to avoid Type I error?

Thanks again!

[Francesca Pesola]

Sorry, one more thing.

if using the Variance-comparison testst (-sdtest-) I found that the variance does not differ across arms, would I be able to use the fully clustered approach described in your paper or would I need additional tests?

[Scott Baldwin]

No, arm should not be missing in the control. Just assign each control subject to his/her own cluster (cluster of 1). If you specify the model, like we describe in the paper, they won't contribute to the estimate of the cluster variance. For the residuals, use the -residual- option of -mixed-.

mixed fY arm Ycentred Centre covs1 cov2 cov3 || clust: arm, nocons vce(robust) residual(ind, by(arm))

You could perform a likelihood ratio test between the model with heteroscedastic residuals (the one right above) and a model with homoscedastic residuals (same as above, just exclude the -residual- option).

mixed fY arm Ycentred Centre covs1 cov2 cov3 || clust: arm, nocons vce(robust) residual(ind, by(arm))
estimates store het

mixed fY arm Ycentred Centre covs1 cov2 cov3 || clust: arm, nocons vce(robust) 
estimates store hom

lrtest het hom

The null is that the heteroscedastic and homoscedastic models provide identical fit.

Best,
Scott

[Francesca Pesola]

Hi Scott,

Thanks again for such a helpful reply. I now code therapist clusters, in the experimental arm, as 0 and 1 and then number each participant in the control group as an independent cluster (2-110).

I have run the analysis as you suggested and have a few more questions I am hoping you can help with:

In your paper, it appears that the fully and partially nested model are specified similarly - with each person in the control arm treated as a cluster - but the fully nested model assumes equal variances across conditions. So is using -residual(ind, by(arm))- specifying a partially nested model?

The LRtest is not significant (p = .21) so I can assume the variances are homogeneous across arms and report the results for the homogenous model, correct? In this case would I be able to conclude that the fully and partially nested model approaches are both correct (or even equivalent)? What should I have done, had the LRtest been significant?

Could -xtmixed- do a similar job?

The variance associated with cluster is very small - could this be ground for dropping it at all?

Random-effects Parameters Estimate Std. Err. [95% Conf. Interval]

clust: Identity
var(arm) 2.56e-15 1.23e-14 2.05e-19 3.21e-11

var(Residual) 249.2436 24.98691 204.7813 303.3596

LR test vs. linear regression: chibar2(01) = 0.00 Prob >= chibar2 = 1.0000

Thanks again,
Francesca

[Scott Baldwin]

The distinction between the fully nested model and partially nested model is separate from the distinction between homo- and heteroscedastic (in our paper, anyway). The former has to do with the random effects and the latter has to do with the residuals. The -residual- option in -mixed- controls the residuals; therefore, just using that option doesn't affect whether the model is fully or partially nested.

With only two clusters, I wouldn't trust the variance component estimate (small or large). I come back to my point in a previous post, you are in a tough spot given the design of the study. No good options.

Best,
Scott

[Francesca Pesola]

Thanks Scott,

I have one more question, in case of a binary outcome I was planning to specify my model using -xtmelogit- as below, is that suffcient?

xtmelogit fCGI_health arm Centre highSF highD highA || clust : arm, nocons or

I am still unsure how mixed (or xtmelogit) would distinguish between a fully and partially nested model, is it by suppressing the overall random intercept while applying a random slope? I'd be grateful if you could point me to any reference or online resource.

Thanks again,
Francesca

[Scott Baldwin]

Looks right to me.

Scott

[Francesca Pesola]

Great - thanks again Scott!