[R] two lme questions
Scott Rifkin
scott.rifkin at yale.edu
Thu Mar 18 23:17:47 CET 2004
This is an update to my previous post today after finding some previous
posts about crossed random effects. Any comments would be much
appreciated.
I have two random factors (Plot and Variety), one fixed factor (Time). I
got rid of one of the fixed factors from my previous post for simplicity.
Each Variety is tested at each of two Times, so I would like to include
Time %in% Variety in my random factors.
I tried the following based on a suggestion from Douglas Bates
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/26424.html
> Data = groupedData(Measurement ~Time|Variety, data =
read.table("mvone.txt",header = TRUE))
> Data$vt = factor(paste(Data$Variety,Data$Time,sep="/"))
> Data$const=factor(rep(1,nrow(Data)))
> Data.lme=lme(Measurement~Time,data=Data,random=list(const=pdBlocked(list(pdIdent(~Plot-1),pdIdent(~Variety-1),pdIdent(~vt-1)))))
I'd like some help interpreting the results.(the ... mean that I deleted
repetitive stuff)
>summary(Data.lme)
...
Block 1: Plot18Ea18La, Plot18Ea49Ea, Plot18Ea59Ea, Plot18Ea67Ea,
Plot18La18Ea, Plot18La21La, Plot18La46La, Plot18La71La, Plot21Ea21La,
...
Plot71La67La, Plot71La71Ea
Formula: ~Plot - 1 | const
Structure: Multiple of an Identity
Plot18Ea18La Plot18Ea49Ea Plot18Ea59Ea Plot18Ea67Ea Plot18La18Ea
StdDev: 0.973655 0.973655 0.973655 0.973655 0.973655
Plot18La21La Plot18La46La Plot18La71La Plot21Ea21La Plot21Ea59Ea
StdDev: 0.973655 0.973655 0.973655 0.973655 0.973655
...
I take it that these StdDevs are the sqrt of the variance component for
Plot. They are the diagonal of the pdBlocked matrix corresponding to
pdIdent(~Plot-1)
same thing for
Block 2: VarietyL23, VarietyL54, VarietyL59, VarietyL71, VarietyL49,
...
Formula: ~Variety - 1 | const
Structure: Multiple of an Identity
VarietyL23 VarietyL54 VarietyL59 VarietyL71 VarietyL49 VarietyL21
StdDev: 0.04296328 0.04296328 0.04296328 0.04296328 0.04296328 0.04296328
...
and
Block 3: vtL18/Early, vtL18/Late, vtL21/Early, vtL21/Late, vtL23/Early,
...
vtL71/Late
Formula: ~vt - 1 | const
Structure: Multiple of an Identity
vtL18/Early vtL18/Late vtL21/Early vtL21/Late vtL23/Early vtL23/Late
StdDev: 0.006307627 0.006307627 0.006307627 0.006307627 0.006307627 0.006307627
...
and
Residual
StdDev: 0.1299986
is just the sigma(e)
Fixed effects: Measurement ~ Time
Value Std.Error DF t-value p-value
(Intercept) 9.762629 0.10228970 190 95.44097 <.0001
TimeLate -0.222656 0.03712913 190 -5.99680 <.0001
I'm still wondering why I only get TimeLate...is it because TimeEarly is
just 0.222656?
when I do:
> intervals(Data.lme)
Approximate 95% confidence intervals
Fixed effects:
lower est. upper
(Intercept) 9.5608597 9.7626290 9.9643984
TimeLate -0.2958942 -0.2226559 -0.1494177
attr(,"label")
[1] "Fixed effects:"
Random Effects:
Level: const
lower est. upper
sd(Plot - 1) 8.436375e-01 0.973655068 1.1237103
sd(Variety - 1) 1.577451e-02 0.042963279 0.1170143
sd(vt - 1) 4.186565e-06 0.006307627 9.5032935
Within-group standard error:
lower est. upper
0.1119908 0.1299986 0.1509019
These are my confidence intervals on the variance components, I assume.
I'm slightly confused on this because in one of the posts I found:
http://finzi.psych.upenn.edu/R/Rhelp02a/archive/21857.html
it seems to be much more complicated to find the variance componenets.
Thanks much,
Scott Rifkin
scott.rifkin at yale.edu
More information about the R-help
mailing list