Hi all,

I have a problem based on a colleague's experiment that I've been asked to analyse, which is more of a general mixed modelling issue rather than specifically an R issue, and I would be extremely grateful for any help that any readers of this list can provide.

An experiment was conducted in which the aim was to compare 3 concentrations of 2 active treatments (i.e. 6 active treatments in total) to a negative control.  Three batches of each of the actives have been tested, and 3 reps tested for each batch.  In contrast, 20 replicates have been taken of the negative control - but, by definition, there is no "batch" for this treatment.

Here is some code to reproduce the experimental design:

Treat<- factor(c(rep("NC", 20), rep("A", 27), rep("B", 27)))
Conc<-factor(c(rep(1, 20), rep(1:3, each=9), rep(1:3, each=9)))
Batch<-factor(c(rep(1, 20), rep( rep(1:3, each=3), 6)))
Treatment<-factor(Treat:Conc)  #specify new treatment variable (so don't attempt to estimate Conc. 2&3 for NC)

I originally planned to analyses these data in a LMM, with Treat*Conc as a 7 level fixed effect (i.e. 3*2 actives + control), and with Treat:Conc:Batch as random.  The following code simulates my response variable assuming this model:

                Resp<-  rep(9, 74) + #simulate intercept
                                c( rep(rnorm(1, 0, sd=2.5), 20)^2, rep(rnorm(18, 0, sd=2.5), each=3)^2) + #simulate treat.conc.batch variance
                                rep(rnorm(74, 0, sd=.2)^2) + #simulate residual variance
                                c(rep(0,20), rep(c(-4, 0,0,-4, 0,0), each= 9)) #simulate fixed effects
                Data<-data.frame(Treatment, Conc, Batch, Resp)

While this code models the data using lme4:
                Mod<-lmer(Resp ~ Treatment + (1|Treatment:Batch), data=Data)

I can now obtain and plot treatment means/CIs using glht in the multcomp package:
library(multcomp)
                Mean.mat<-diag(rep(1,7))
                                Mean.mat[,1]<-rep(1,7)
                                rownames(Mean.mat)<-levels(Data$Treatment)
                Est.means<-glht(Mod, Mean.mat)
plot(Est.means)

Hopefully from the above plot you can see what my issue is.  The negative control, which I want to compare everything against, has by far the least precision around its estimate, despite the data for the control hardly varying at all.  This happens because the greatest source of variability in the model (by far) is the variability between batches, but different batches of the negative control don't exist.  As such, I'm not sure that this is a fair way to model the data, because the negative control is unfairly penalised by the variability between the batches of the other treatments.

I imagine that this kind of problem isn't particularly uncommon, but it's the first time I've had to deal with something like this myself.   The only potential solution I've come up with so far is to scrap the negative control from the model, and simply subtract the negative control's mean "count" from all other values (either by specifying this mean as an offset or by subtracting it from all data-points).  But this will probably give "anti-conservative" results, as it would assume the mean for the negative control doesn't vary.

I would be extremely grateful if anyone would care to share their thoughts on possible solutions to this problem - and whether anyone has dealt with this kind of issue before.  I feel that I may well be missing something obvious - but can't see at the moment how else to get around it!

Many thanks for any help you can provide.

Cheers,

Nick





	[[alternative HTML version deleted]]


