[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