[R] lme syntax for P&B examples

Henric Nilsson henric.nilsson at statisticon.se
Thu Feb 9 14:25:17 CET 2006


Paul Cossens said the following on 2006-02-09 06:21:
> Hi Harold,
> 
> 
> Thanks for your reply. I had already looked at all the reading material
> you suggested but updated to the latest Matrix 
> as recommneded then spent all day trying to figure out what is
> happening.  
> 
> I worked through the problems and give my workings below that others may
> find useful. 
> (My notation is to use lme> to show lme commands and lmer> to show lmer
> commands. 
> I worked on two sessions in parallel. My comments are preceded by double
> hashes '##' and
> questions '##??'. I haven't included the datasets.)   
> 
> I have a couple of comments and outstanding issues:
> 
> 1. In the Pixel data set and formulas I think the formulas are printed
> incorrectly in the 
> book as some use 'I(day^2)' while others use just 'day^2'. I have used
> 'I(day^2)'. I'm not sure why the I() function is used. In the fm4Pixel
> example below the answers don't match up exactly but are close.
> 
> The lme example is 
> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> which I have converted to lmer:
> fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data = Pixel)
> 
> The t-values for Side are close (sse below) but different enough to
> wonder if I am still doing something wrong? 

It's wrong. Try

fm6Pixel <- lmer(pixel ~ day + I(day^2) + Side + (day|Dog) + 
(1|Side:Dog), data = Pixel)

> 2. To me the specification description in the R-News article is
> confusing as it seems 
> to suggest that nesting does not need to be completely specified if the
> groupings and nestings are clear in data set. 
> 
> Prof Bates article in R news vol 5/1 P 30  states "It happens in this
> case that the grouping factors 'id' and 'sch' are not nested but if they
> were nested there would be no change in the model specification"
> 
> If the lme formula is 
> fm1Oxide<-lme(Thickness~1,Oxide)
> 
> I have found the formula lmer parlance should be:
> 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)'
> not 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)'
> as the article reads to me. 
> 
> In other words you always need to explicitly specify nesting levels.
> 
> 3. I still can't figue out how to replicate the lme formula
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> i.e 
> formula(fm2Oxide)
> Thickness ~ 1
> 
> If I simply drop the Lot:Wafer term as in 'fm2Oxide<-lmer(Thickness~
> (1|Lot),data=Oxide)'
> I get the error
> 
> 'Error in x[[3]] : object is not subsettable'
> 
> what's the solution?

fm3Oxide <- lmer(Thickness ~ 1 + (1|Lot), data = Oxide)

This seems to be a bug in `lmer', since one doesn't have to specifiy the 
intercept explicitly in e.g.

lmer(Thickness ~ (1|Lot) + (1|Wafer), data = Oxide)


HTH,
Henric



> 
> I'd be interested to read you article for further insights.
> 
> Thanks
> 
> Paul
> 
> 
> 
> #############################################################
> #Oxide
> # P&B(2000) p167-170
> 
> #NLME lme example
> 
> lme>data(Oxide)
> lme>formula(Oxide)
> lme>Thickness ~ 1 | Lot/Wafer
> lme>fm1Oxide<-lme(Thickness~1,Oxide)
> lme> fm1Oxide
> Linear mixed-effects model fit by REML
>   Data: Oxide 
>   Log-restricted-likelihood: -227.0110
>   Fixed: Thickness ~ 1 
> (Intercept) 
>    2000.153 
> 
> Random effects:
>  Formula: ~1 | Lot
>         (Intercept)
> StdDev:    11.39768
> 
>  Formula: ~1 | Wafer %in% Lot
>         (Intercept) Residual
> StdDev:    5.988802 3.545341
> 
> Number of Observations: 72
> Number of Groups: 
>            Lot Wafer %in% Lot 
>              8             24 
> 
> lme> intervals(fm1Oxide, which = "var-cov")
> Approximate 95% confidence intervals
> 
>  Random Effects:
>   Level: Lot 
>                   lower     est.    upper
> sd((Intercept)) 6.38881 11.39768 20.33355
>   Level: Wafer 
>                    lower     est.    upper
> sd((Intercept)) 4.063919 5.988802 8.825408
> 
>  Within-group standard error:
>    lower     est.    upper 
> 2.902491 3.545341 4.330572 
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> lme> fm2Oxide
> Linear mixed-effects model fit by REML
>   Data: Oxide 
>   Log-restricted-likelihood: -245.5658
>   Fixed: Thickness ~ 1 
> (Intercept) 
>    2000.153 
> 
> Random effects:
>  Formula: ~1 | Lot
>         (Intercept) Residual
> StdDev:    11.78447 6.282416
> 
> Number of Observations: 72
> Number of Groups: 8
> 
> lme>intervals(fm2Oxide, which = "var-cov")
> Approximate 95% confidence intervals
> 
>  Random Effects:
>   Level: Lot 
>                    lower     est.    upper
> sd((Intercept)) 6.864617 11.78447 20.23035
> 
>  Within-group standard error:
>    lower     est.    upper 
> 5.283116 6.282416 7.470733 
> 
> lme> coef(fm1Oxide, level=1)
>   (Intercept)
> 1    1996.689
> 2    1988.931
> 3    2001.022
> 4    1995.682
> 5    2013.616
> 6    2019.561
> 7    1991.954
> 8    1993.767
> coef(fm1Oxide, level=1)
>   (Intercept)
> 1    1996.689
> 2    1988.931
> 3    2001.022
> 4    1995.682
> 5    2013.616
> 6    2019.561
> 7    1991.954
> 8    1993.767
>> coef(fm1Oxide, level=2)
>     (Intercept)
> 1/1    2003.235
> 1/2    1984.730
> 1/3    2001.146
> 2/1    1989.590
> 2/2    1988.097
> 2/3    1986.008
> 3/1    2002.495
> 3/2    2000.405
> 3/3    2000.405
> 4/1    1995.668
> 4/2    1998.951
> 4/3    1991.191
> 5/1    2009.184
> 5/2    2016.646
> 5/3    2018.735
> 6/1    2031.296
> 6/2    2021.745
> 6/3    2011.000
> 7/1    1990.204
> 7/2    1991.398
> 7/3    1991.995
> 8/1    1993.677
> 8/2    1995.170
> 8/3    1990.693
> 
> 
> ######## the is the lmer example using Matrix 0.995-5
> 
> 
> lmer>Oxide<-read.csv("Oxide.csv",header=TRUE);
> lmer>Oxide$Source<-as.factor(Oxide$Source)
> lmer>Oxide$Lot<-as.factor(Oxide$Lot)
> lmer>Oxide$Wafer<-as.factor(Oxide$Wafer)
> Lmer>Oxide$Site<-as.factor(Oxide$Site)
> 
> 
> ## Bates article in R news vol5/1 says specifying nesting explicity
> isn't necessary: 
> ## P 30 "It happens in this case that the grouping factors 'id' and
> 'sch' are not 
> ## nested but if they were nested there would be no change in the model
> specification"
> ## Following this one would expect that the following statement would
> automatically 
> ## detect nesting
> 
> lmer>fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)
> 
> 
> lmer> fm1Oxide
> Linear mixed-effects model fit by REML
> Formula: Thickness ~ (1 | Lot) + (1 | Wafer) 
>    Data: Oxide 
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  496.6093 503.4393 -245.3046   495.3528     490.6093
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  Lot      (Intercept) 138.9981 11.7897 
>  Wafer    (Intercept)   1.4930  1.2219 
>  Residual              38.3490  6.1927 
> # of obs: 72, groups: Lot, 8; Wafer, 3
> 
> Fixed effects:
>              Estimate Std. Error t value
> (Intercept) 2000.1528     4.2901  466.22
> 
> ## The lme vs lmer std.devs for Lot are 11.39768 : 11.7987  
> ## The lme vs lmer std.devs for Wafer are 5.988802 : 1.2219 
> ## If my lmer specifcation is correct then the Wafer std.dev seems too
> big.
> ## also the levels aren't specifed correctly as the following ranef()
> function
> ## shows
> 
> lmer> ranef(fm1Oxide)
> An object of class "lmer.ranef"
> [[1]]
>   (Intercept)
> 1  -3.7058415
> 2 -12.0069264
> 3   0.9298293
> 4  -4.7839045
> 5  14.4056166
> 6  20.7661881
> 7  -8.7727375
> 8  -6.8322241
> 
> [[2]]
>   (Intercept)
> 1   0.9526407
> 2  -0.2750582
> 3  -0.6775825
> 
> ###There is no nesting of wafers within lots
> ###Note however that the following appears to work the same as in the
> P&B text
>  
> lmer> fm3Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)
> lmer> fm3Oxide
> Linear mixed-effects model fit by REML
> Formula: Thickness ~ (1 | Lot) + (1 | Lot:Wafer) 
>    Data: Oxide 
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  460.0221 466.8521 -227.0110   458.7378     454.0221
> Random effects:
>  Groups    Name        Variance Std.Dev.
>  Lot:Wafer (Intercept)  35.866   5.9888 
>  Lot       (Intercept) 129.853  11.3953 
>  Residual               12.570   3.5454 
> # of obs: 72, groups: Lot:Wafer, 24; Lot, 8
> 
> Fixed effects:
>              Estimate Std. Error t value
> (Intercept) 2000.1528     4.2309  472.75
> 
> ## the following gives the coefficients for levels however the number of
> levels is 
> ## in reverse order to lme
> ## Also note that the order of levels seems to have changed from the
> lmer model 
> ## fm1Oxide where Lot coefficnets are in [[1]] and Wafer [[2]] 
> lmer>ranef(fm3Oxide)
> An object of class "lmer.ranef"
> [[1]]
>      (Intercept)
> 1:1   6.54582998
> 1:2 -11.95898390
> 1:3   4.45657680
> 2:1   0.65819690
> 2:2  -0.83412680
> 2:3  -2.92337998
> 3:1   1.47284009
> 3:2  -0.61641309
> 3:3  -0.61641309
> 4:1  -0.01366505
> 4:2   3.26944709
> 4:3  -4.49063615
> 5:1  -4.43133801
> 5:2   3.03028049
> 5:3   5.11953367
> 6:1  11.73559478
> 6:2   2.18472310
> 6:3  -8.56000754
> 7:1  -1.74970878
> 7:2  -0.55584982
> 7:3   0.04107966
> 8:1  -0.09041889
> 8:2   1.40190481
> 8:3  -3.07506629
> 
> [[2]]
>   (Intercept)
> 1   -3.463334
> 2  -11.221203
> 3    0.868982
> 4   -4.470850
> 5   13.462925
> 6   19.407266
> 7   -8.198657
> 8   -6.385129
> 
> ##?? given that the fm3Oxide works one would expect that dropping the 
> ##?? (1|Lot:Wafer) term would work however it give the following error
> 
> lmer>fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide)
> Error in x[[3]] : object is not subsettable
> 
> ##?? the question is how to specify the lme model
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> ##?? in lme syntax
> 
> 
> ######################################################################
> 
> #Pixel
> # P&B(2000) p40-45
> 
> ## the lme example is as follows
> 
> lme>data(Pixel)
> 
> ## firstly the book may have an error on P 42 line 1
> lme_wrong> fm1Pixel<-lme(pixel~day+day^2,data=Pixel,random = list(Dog= ~
> day ,Side= ~ 1))
> ###which gives: 
> lme_wrong> intervals(fm1Pixel)
> Approximate 95% confidence intervals
> 
>  Fixed effects:
>                   lower       est.       upper
> (Intercept) 1071.415167 1093.21538 1115.015590
> day           -1.126053   -0.14867    0.828713
> attr(,"label")
> [1] "Fixed effects:"
> 
>  Random Effects:
>   Level: Dog 
>                           lower       est.      upper
> sd((Intercept))      17.3485293 31.4936604 57.1720308
> sd(day)               0.3586459  1.0719672  3.2040342
> cor((Intercept),day) -0.9822311 -0.7863546  0.2294876
>   Level: Side 
>                    lower     est.    upper
> sd((Intercept)) 8.519543 15.08995 26.72755
> 
>  Within-group standard error:
>    lower     est.    upper 
> 12.35975 14.53391 17.09052 
> 
> ## the book should probably give I(day^2) instead of day^2 on p42 line 1
> where 
> ## as this gives the answer in the book
> lme> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> lme> intervals(fm1Pixel)
> Approximate 95% confidence intervals
> 
>  Fixed effects:
>                    lower         est.        upper
> (Intercept) 1053.0968388 1073.3391382 1093.5814376
> day            4.3796925    6.1295971    7.8795016
> I(day^2)      -0.4349038   -0.3673503   -0.2997967
> attr(,"label")
> [1] "Fixed effects:"
> 
>  Random Effects:
>   Level: Dog 
>                           lower       est.      upper
> sd((Intercept))      15.9284469 28.3699038 50.5291851
> sd(day)               1.0812133  1.8437505  3.1440751
> cor((Intercept),day) -0.8944371 -0.5547222  0.1909581
>   Level: Side 
>                    lower     est.    upper
> sd((Intercept)) 10.41726 16.82431 27.17195
> 
>  Within-group standard error:
>     lower      est.     upper 
>  7.634522  8.989606 10.585209 
> 
> lme>VarCorr(fm1Pixel) 
>             Variance       StdDev    Corr  
> Dog =       pdLogChol(day)                 
> (Intercept) 804.851443     28.369904 (Intr)
> day           3.399416      1.843750 -0.555
> Side =      pdLogChol(1)                   
> (Intercept) 283.057248     16.824305       
> Residual     80.813009      8.989606
> 
> lme>  summary(fm1Pixel)
> Linear mixed-effects model fit by REML
>  Data: Pixel 
>        AIC      BIC    logLik
>   841.2102 861.9712 -412.6051
> 
> Random effects:
>  Formula: ~day | Dog
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev    Corr  
> (Intercept) 28.369904 (Intr)
> day          1.843750 -0.555
> 
>  Formula: ~1 | Side %in% Dog
>         (Intercept) Residual
> StdDev:    16.82431 8.989606
> 
> Fixed effects: pixel ~ day + I(day^2) 
>                 Value Std.Error DF   t-value p-value
> (Intercept) 1073.3391 10.171686 80 105.52225       0
> day            6.1296  0.879321 80   6.97083       0
> I(day^2)      -0.3674  0.033945 80 -10.82179       0
>  Correlation: 
>          (Intr) day   
> day      -0.517       
> I(day^2)  0.186 -0.668
> 
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max 
> -2.82905723 -0.44918107  0.02554930  0.55721629  2.75196509 
> 
> Number of Observations: 102
> Number of Groups: 
>           Dog Side %in% Dog 
>            10            20  
> 
> ## Note thate the formula in the book is on P44 sect 1.5.1 summary table
> is
> ##  Fixed effects: pixel ~ day + day^2 
> ## compared to above:
> ## Fixed effects: pixel ~ day + I(day^2)
> ## but the anova on page 45 is the same
> 
> lme>fm2Pixel <- update(fm1Pixel,random = ~day | Dog)
> lme>anova(fm1Pixel,fm2Pixel)
>          Model df      AIC      BIC    logLik   Test L.Ratio p-value
> fm1Pixel     1  8 841.2102 861.9712 -412.6051                       
> fm2Pixel     2  7 884.5196 902.6854 -435.2598 1 vs 2 45.3094  <.0001
> 
> 
> lme> fm3Pixel <- update(fm1Pixel,random = ~1 | Dog/Side)
> lme > anova(fm1Pixel,fm3Pixel)
>          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> fm1Pixel     1  8 841.2102 861.9712 -412.6051                        
> fm3Pixel     2  6 876.8390 892.4098 -432.4195 1 vs 2 39.62885  <.0001
> 
> 
> ##again there following seems to be a mistake in the book
> 
> lme> fm4Pixel <- update(fm1Pixel,pixel ~ day + day^2 + Side)
> lme> summary(fm4Pixel)
> Linear mixed-effects model fit by REML
>  Data: Pixel 
>       AIC     BIC    logLik
>   895.071 915.832 -439.5355
> 
> Random effects:
>  Formula: ~day | Dog
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev    Corr  
> (Intercept) 31.520129 (Intr)
> day          1.073342 -0.786
> 
>  Formula: ~1 | Side %in% Dog
>         (Intercept) Residual
> StdDev:    15.01697 14.50742
> 
> Fixed effects: pixel ~ day + Side 
>                 Value Std.Error DF  t-value p-value
> (Intercept) 1097.5272 11.562590 81 94.92053  0.0000
> day           -0.1496  0.491218 81 -0.30451  0.7615
> SideR         -8.6098  7.379984  9 -1.16664  0.2733
>  Correlation: 
>       (Intr) day   
> day   -0.633       
> SideR -0.319  0.000
> 
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max 
> -3.73906417 -0.38367706  0.04758941  0.39690056  2.23720545 
> 
> Number of Observations: 102
> Number of Groups: 
>           Dog Side %in% Dog 
>            10            20 
> 
> ###the book should probably give I(day^2) instead of day^2
> ### as the fixed effect coefficient names in the summary table 
> ##  differ from the formula given
> ##  the following seems a partial correction   
> lme> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> lme>summary(fm5Pixel)
> 
> Linear mixed-effects model fit by REML
>  Data: Pixel 
>        AIC      BIC    logLik
>   835.8546 859.1193 -408.9273
> 
> Random effects:
>  Formula: ~day | Dog
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev    Corr  
> (Intercept) 28.463605 (Intr)
> day          1.843823 -0.553
> 
>  Formula: ~1 | Side %in% Dog
>         (Intercept) Residual
> StdDev:     16.5072 8.983614
> 
> Fixed effects: pixel ~ day + I(day^2) + Side 
>                 Value Std.Error DF   t-value p-value
> (Intercept) 1077.9484 10.862705 80  99.23388  0.0000
> day            6.1296  0.879023 80   6.97323  0.0000
> I(day^2)      -0.3674  0.033923 80 -10.82914  0.0000
> SideR         -9.2175  7.626768  9  -1.20858  0.2576
>  Correlation: 
>          (Intr) day    I(d^2)
> day      -0.484              
> I(day^2)  0.174 -0.667       
> SideR    -0.351  0.000  0.000
> 
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max 
> -2.80982455 -0.47133415  0.02610263  0.54115378  2.77470104 
> 
> Number of Observations: 102
> Number of Groups: 
>           Dog Side %in% Dog 
>            10            20 
> 
> 
> ## Note however that the Fixed effects Values above differ from the
> values in the book  
> 
> ####### the is the lmer example for Pixel using Matrix 0.995-5
> 
> ##the lmer example for the Pixel data is exported and reimported as a
> csv file
> lmer>Pixel<-read.csv("Pixel.csv",header=TRUE);
> lmer>Pixel$Side<-as.factor(Pixel$Side)
> lmer>Pixel$Dog<-as.factor(Pixel$Dog)
> lmer>Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
>  
> lmer>fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Dog:Side),
> data = Pixel)
> lmer>fm1Pixel
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side) 
>    Data: Pixel 
>       AIC     BIC    logLik MLdeviance REMLdeviance
>  839.2102 857.585 -412.6051   827.3298     825.2102
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr   
>  Dog:Side (Intercept) 283.0567 16.8243         
>  Dog      (Intercept) 804.8460 28.3698         
>           day           3.3994  1.8438  -0.555 
>  Residual              80.8130  8.9896         
> # of obs: 102, groups: Dog:Side, 20; Dog, 10
> 
> Fixed effects:
>                Estimate  Std. Error t value
> (Intercept) 1073.339126   10.171658 105.523
> day            6.129599    0.879321   6.971
> I(day^2)      -0.367350    0.033945 -10.822
> 
> Correlation of Fixed Effects:
>          (Intr) day   
> day      -0.517       
> I(day^2)  0.186 -0.668
> 
> ## 
> 
> lme> VarCorr(fm1Pixel)
> $"Dog:Side"
> 1 x 1 Matrix of class "dpoMatrix"
>             (Intercept)
> (Intercept)    283.0567
> 
> $Dog
> 2 x 2 Matrix of class "dpoMatrix"
>             (Intercept)        day
> (Intercept)   804.84605 -29.015418
> day           -29.01542   3.399415
> 
> attr(,"sc")
> [1] 8.989606
> 
> lmer> fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel)
> lmer> fm2Pixel
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + (day | Dog) 
>    Data: Pixel 
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  882.5196 898.2694 -435.2598   873.5964     870.5196
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr   
>  Dog      (Intercept) 892.7720 29.8793         
>           day           3.0772  1.7542  -0.488 
>  Residual             197.5767 14.0562         
> # of obs: 102, groups: Dog, 10
> 
> Fixed effects:
>               Estimate Std. Error t value
> (Intercept) 1072.92725   10.44452 102.726
> day            6.08910    1.14699   5.309
> I(day^2)      -0.35677    0.05221  -6.833
> 
> Correlation of Fixed Effects:
>          (Intr) day   
> day      -0.541       
> I(day^2)  0.286 -0.799
> 
> lmer> anova(fm1Pixel,fm2Pixel)
> 
> Data: Pixel
> Models:
> fm2Pixel: pixel ~ day + I(day^2) + (day | Dog)
> fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
>          Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)    
> fm2Pixel  6  882.52  898.27 -435.26                             
> fm1Pixel  7  839.21  857.59 -412.61 45.309      1  1.682e-11 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> ## Some of the statistics here are slightly different from above,
> notably the Df 
> ## but I guess the result is the same
> 
> lmer> fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data =
> Pixel)
> lmer> anova(fm1Pixel,fm3Pixel)
> 
> Data: Pixel
> Models:
> fm3Pixel: pixel ~ day + I(day^2) + (1 | Dog:Side)
> fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
>          Df     AIC     BIC  logLik Chisq Chi Df Pr(>Chisq)    
> fm3Pixel  4  877.88  888.38 -434.94                            
> fm1Pixel  7  839.21  857.59 -412.61 44.67      3  1.087e-09 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> 
> ## Some of the statistics are slightly different again
> 
> lmer>fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data =
> Pixel)
> lmer>fm4Pixel
> 
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + Side + (day | Dog) 
>    Data: Pixel 
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  876.8204 895.1952 -431.4102   869.6765     862.8204
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr   
>  Dog      (Intercept) 896.1278 29.9354         
>           day           3.0972  1.7599  -0.490 
>  Residual             190.9227 13.8175         
> # of obs: 102, groups: Dog, 10
> 
> Fixed effects:
>                Estimate  Std. Error t value
> (Intercept) 1075.649999   10.521427 102.234
> day            6.091506    1.133983   5.372
> I(day^2)      -0.357334    0.051369  -6.956
> SideR         -5.401961    2.736268  -1.974
> 
> Correlation of Fixed Effects:
>          (Intr) day    I(d^2)
> day      -0.535              
> I(day^2)  0.279 -0.795       
> SideR    -0.130  0.000  0.000
> 
> ##?? Fixed effects estimates are sligtly different 
> ##?? As df and p-values are missing I assume that it can be concluded
> that as 
> ##?? t-value is less than 1.96 that 'Side' is not sigificant. 
> ##?? In the lme example for fm4Pixel the t-value for Side is -1.21
> ##?? have i specified fm4Pixel correctly?
> 
> -----Original Message-----
> From: Doran, Harold [mailto:HDoran at air.org] 
> Sent: Thursday, 9 February 2006 02:01
> To: Paul Cossens; r-help at stat.math.ethz.ch
> Subject: RE: [R] lme syntax for P&B examples
> 
> 
> Paul:
> 
> It is a little difficult to understand what you are trying to translate
> since you do not show what the model would look like using lme. If you
> show lme, then it is easy to translate into lmer syntax.
> 
> A few thoughts, first, use lmer in the Matrix package and not in lme4.
> Second, see the Bates article in R news at the link below for dealing
> with nesting structures. Last, a colleague and I have a paper in press
> showing how to fit models using lme which we submitted a year or so ago.
> Since lme has evolved to lmer, we created an appendix that translates
> all of our lme models to the lmer syntax so readers can see
> equivalences. I am happy to send this to you (or others) upon request.
> 
> http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf
> 
> Harold
> 
> 
>  
> 
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Paul Cossens
> Sent: Wednesday, February 08, 2006 12:08 AM
> To: r-help at stat.math.ethz.ch
> Subject: [R] lme syntax for P&B examples
> 
> Hi helpeRs,
>  
> I've been working through some examples in Pinhiero & Bates( 2000)
> trying to understand how to translate to the new Lme4 syntax but without
> much luck. Below is what I think I should do, but either the answers
> don't come out the same or I get errors. 
> In the Oxide problems I'm particularly interested in obtaining the
> levels coeficients but this options no longer seems to be available in
> lme4. How can levels infor be obtained in lme4?
>  
> If someone can recreate the examples below in lme4 syntax so I can
> follow what is happening in the text I'd be grateful. 
>  
> Cheers
>  
> Paul Cossens
>  
>  
> #Pixel
> # P&B(2000) p40-45
>  
> Pixel<-read.csv("Pixel.csv",header=TRUE);
> Pixel$Side<-as.factor(Pixel$Side)
> Pixel$Dog<-as.factor(Pixel$Dog)
>  
> (fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Side), data =
> Pixel))
> (fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel))
> (fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data = Pixel))
> or should I do it this way? Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
> (fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|DS), data = Pixel))
>  
> (fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side , data = Pixel))
>  
> 
> #Oxide
> # P&B(2000) p167-170
>  
> Oxide<-read.csv("Oxide.csv",header=TRUE);
> Oxide$Source<-as.factor(Oxide$Source)
> Oxide$Lot<-as.factor(Oxide$Lot)
> Oxide$Wafer<-as.factor(Oxide$Wafer)
> Oxide$Site<-as.factor(Oxide$Site)
> fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide) )
> (fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide) )
> coef(fm1Oxide)
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html




More information about the R-help mailing list