On Mon, Apr 12, 2010 at 8:00 AM, Kevin E. Thorpe
<kevin.thorpe at utoronto.ca> wrote:
Kevin E. Thorpe <kevin.thorpe at ...> writes:
My data come from a crossover trial and are balanced.
'data.frame': 96 obs. of 4 variables:
$ Subject : int 1 2 3 5 6 7 10 11 12 13 ...
$ Treatment: Factor w/ 2 levels "Barley","Oat": 1 1 1 1 1 1 1
1 1 1 ... $ Dose : int 8 8 8 8 8 8 8 8 8 8 ...
$ iAUC : num 110 256 129 207 244 ...
clip>
Shouldn't you make Subject into a factor?
Ken
It would make the plot a little bit prettier but I don't think it
matters in this case because variable that appears as a grouping
variable (i.e. on the right of the | ) is automatically treated
as a factor? I think?
Since it is really a crossover trial, it would seem reasonable in
principle to have the (Treatment|Subject) random effect in there
as well. I'm not sure what to do about the -1 correlation: it
seems the choices (not necessarily in order) are (1) throw up
your hands and say there's not enough data to estimate
independently; (2) try WinBUGS, possibly with slightly
informative priors; (3) try using lme4a to create profiles of the
parameters and see if you can figure out what's happening.
Let's see. I wish (1) was an option. (2) would be promising if my
knowledge of BUGS and Bayesian methods filled more than a thimble.
Thanks to Jarrod for his suggestion in response to this. I'll take
a look at that too. Option (3) is probably worth a go too.
Aside from the fact that the Dose variable are the actual doses and
not categories, and we all know not to categorize continuous
variables, what are your thoughts on treating Dose as a factor
(since it seems to behave)?
Thanks all for taking the time to provide your suggestions.
Kevin
Regarding lme4a: how do I obtain it? I guess that comes down to,
what is the repository to give to install.packages()? Does it
require a different Matrix package than the one I have, which is,
0.999375-33 and if so, how do I not break my current lme4/Matrix
combination?
The repository is http://R-forge.R-project.org but be aware that
lme4a is under active development and not guaranteed against
breakage. It would be inadvisable to rely on functions and classes
in that package to persist.
By the way, the problems with these data get stranger. In a
different outcome from the same trial the following results from a
model fitting attempt.
lmer(iAUC~Treatment+Dose+(Treatment|Subject)+(Dose|Subject),data=insu
lin)
Linear mixed model fit by REML
Formula: iAUC ~ Treatment + Dose + (Treatment | Subject) + (Dose |
Subject) Data: insulin AIC BIC logLik deviance REMLdev
1956 1982 -968 1983 1936
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 4.2678e-02 2.0659e-01
TreatmentOat 1.6115e+07 4.0144e+03 0.000
Subject (Intercept) 3.0430e+08 1.7444e+04
Dose 1.5173e+06 1.2318e+03 -1.000
And did you notice that you have an Intercept term by Subject in
there twice? It is not surprising that the parameter estimates are
unstable. You will need to rethink the model.
Thanks Doug. Any suggestions? Note that the instability is present
without the random effects for treatment in the model too.
Residual 3.1907e+07 5.6486e+03
Number of obs: 96, groups: Subject, 12
Fixed effects:
Estimate Std. Error t value
(Intercept) 40142.4 5146.7 7.800
TreatmentOat 1340.3 1634.7 0.820
Dose -2675.1 405.5 -6.597
Correlation of Fixed Effects:
(Intr) TrtmnO
TreatmentOt -0.079
Dose -0.922 0.000
As you can see, I get a 0 correlation within one set of random
effects and -1.0 in the other. Also, the fact the fixed effects
estimates are huge makes me suspicious.
Note that if I drop the treatment portions and fit the Dose model to
only one treatment, the correlation is again -1.0 and the fixed
effects are similar.
--
Kevin E. Thorpe
Biostatistician/Trialist, Knowledge Translation Program Assistant
Professor, Dalla Lana School of Public Health University of Toronto
email: kevin.thorpe at utoronto.ca Tel: 416.864.5776 Fax: 416.864.3016