
Hi,

I have been struggling with the variance structure of random effects using lme for a while now and wondered whether anyone could tell me
whether it is possible to fit the model below?

I am fitting a multi-level meta-analysis-style model:

Let k be the clinical indication k=1,...16 (where k=1 if patients in the trial are alcohol dependent, k=2 if patients have binge eating disorder, k=3 if the patients are bipolar disorder, k=4 if the patients have borderline personality disorder,..etc)
Let j denote the trial j=1,....,NS_k (NS_k=number of studies of indication k)
All trials are comparing the same two drugs.

The model I need to fit is
Risk difference_jk = delta + v_k + u_jk
where v_k~Normal(0,gamma^2)
and
u_jk ~ Normal(0, tausq_k)

I am struggling to find a way to allow variance of u_jk (i.e. tausq_k) to differ for each indication. Is this possible?

I have fitted the model that assumes the tausq_k are equal for each indication (ie. tausq_k= tausq) using the data and code below:

Any help or suggestions would be very gratefully received.
many thanks,
Sarah


#Example data
#The risk difference of getting a headache comparing 2 drugs (one risk difference per trial)
Riskdifference=c(-0.08, -0.01, 0.17, -0.02, -0.02, -0.04, 0.00, -0.03, 0.03, 0.03, -0.08, -0.03, -0.14, -0.07, 0.01, 0.04, -0.05, -0.05, -0.08, 0.17, -0.09, 0.04, 0.04, 0.02, 0.11, 0.01, 0.02, -0.11, 0.11, -0.10, 0.00, -0.13, 0.04, 0.07, -0.01)

#The standard error of the risk difference (one per trial)
Stderror=c(0.05, 0.12, 0.12, 0.03, 0.05, 0.08, 0.07, 0.08, 0.07, 0.07, 0.09, 0.03, 0.12, 0.08, 0.08, 0.04, 0.08, 0.08, 0.11, 0.10, 0.08, 0.04, 0.04, 0.05, 0.14, 0.03, 0.03, 0.10, 0.15, 0.13, 0.02, 0.11, 0.07, 0.08, 0.15)

#there are 35 distinct trials
trialno=as.factor(c(seq(1:35)))

#the trials included patients with 16 different clinical indications e.g. epilepsy, migraine, bulimia.
indicationno=as.factor(c(1,2,2,2,2,3,4,5,5,6,7,7,8,8,8,8,8,8,8,8,8,8,9,10,11,11,11,12,13,13,14,14,14,15,16))

#My model assuming fixed tausq
fit = lme(Riskdifference ~ 1, random =~1|indicationno/trialno, method="REML")
summary(fit)

#My model allowing tausq to vary across indications (with error message)
#creating dummary variables indicating which trials are in which indication groups.
trialind1=as.factor(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind2=as.factor(c(0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind3=as.factor(c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind4=as.factor(c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind5=as.factor(c(0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind6=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind7=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind8=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind9=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind10=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
trialind11=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0))
trialind12=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0))
trialind13=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0))
trialind14=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0))
trialind15=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0))
trialind16=as.factor(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1))

const = factor(rep(1, length(trialno)))
groupdat=groupedData(Riskdifference~1|const)
fit = lme(Riskdifference ~ 1, data=groupdat, random = pdBlocked(list(pdIdent(form=~indicationno-1), pdIdent(form=~trialind1-1), pdIdent(form=~trialind2-1),pdIdent(form=~trialind3-1),pdIdent(form=~trialind4-1),pdIdent(form=~trialind5-1),pdIdent(form=~trialind6-1),pdIdent(form=~trialind7-1),pdIdent(form=~trialind8-1),pdIdent(form=~trialind9-1),pdIdent(form=~trialind10-1),pdIdent(form=~trialind11-1),pdIdent(form=~trialind12-1),pdIdent(form=~trialind13-1),pdIdent(form=~trialind14-1),pdIdent(form=~trialind15-1),pdIdent(form=~trialind16-1))))

Error message.....
'Warning message:
In lme.formula(Riskdifference ~ 1, data = groupdat, random = pdBlocked(list(pdIdent(form = ~indicationno -  :
  Fewer observations than random effects in all level 1 groups'

(I will allow for the meta-analytic weights and correct for the sampling errors in s-plus using lmeControl).

	[[alternative HTML version deleted]]

