[R-sig-ME] Random vs. fixed effects

[Robert A LaBudde]

As an addendum, back in the "Old Days", before there were good mixed 
model programs and all you could fit were fixed effect models, you 
had to estimate random effects by fitting the fixed effect model and 
then taking mean squares of the random effect levels as the estimate 
of the variance.

As an example, consider a simple complete randomized block design:

 > require('nlme')
 > set.seed(1)
 > nr<- 2  #reps per case
 > nB<- 4  #blocks
 > nT<- 2  #treatments
 > n<- nT*nB*nr  #total data
 > T<- rep(c(-1,1),nB*nr)
 > B<- factor(c(rep(1,nr*nT),rep(2,nr*nT),rep(3,nr*nT),rep(4,nr*nT)))
 > rB<- rnorm(nB,0,4)[B]
 > e<- rnorm(n,0,1)  #residuals
 > y<- T + rB + e
 > fit4m<- lme(y ~ T, random=~1|B)
 > summary(fit4m)
Linear mixed-effects model fit by REML
  Data: NULL
        AIC      BIC    logLik
   63.22191 65.77814 -27.61096

Random effects:
  Formula: ~1 | B
         (Intercept)  Residual
StdDev:    4.767061 0.8488003

Fixed effects: y ~ T
                 Value Std.Error DF  t-value p-value
(Intercept) 0.5351939 2.3929577 11 0.223654  0.8271
T           0.6916997 0.2122001 11 3.259658  0.0076
  Correlation:
   (Intr)
T 0

Standardized Within-Group Residuals:
         Min          Q1         Med          Q3         Max
-1.76880069 -0.62080393 -0.02900206  0.78562378  1.43929227

Number of Observations: 16
Number of Groups: 4
 > fit4f<- lm(y ~ T + B)
 > summary(fit4f)

Call:
lm(formula = y ~ T + B)

Residuals:
      Min       1Q   Median       3Q      Max
-1.46741 -0.50294 -0.03867  0.66197  1.25562

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  -2.3221     0.4244  -5.472 0.000194 ***
T             0.6917     0.2122   3.260 0.007604 **
B2            3.5997     0.6002   5.998 8.96e-05 ***
B3           -1.4594     0.6002  -2.432 0.033319 *
B4            9.2889     0.6002  15.477 8.20e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8488 on 11 degrees of freedom
Multiple R-squared: 0.9727,     Adjusted R-squared: 0.9628
F-statistic: 98.03 on 4 and 11 DF,  p-value: 1.587e-08

 > d<- coef(fit4f)[3:5]  #block coeffs
 > cB<- c(-0.25*d[1]-0.25*d[2]-0.25*d[3],
+   0.75*d[1]-0.25*d[2]-0.25*d[3],
+   -0.25*d[1]+0.75*d[2]-0.25*d[3],
+   -0.25*d[1]-0.25*d[2]+0.75*d[3]) #convert to mean = 0 basis
 > sd(cB)  #estimate random effect
[1] 4.785915

Note that the two approaches give essentially identical estimates.


At 02:11 PM 4/23/2010, Ben Bolker wrote:
>Here's my question for the group:   Given that it is a reasonable 
>*philosophical* position to say 'treat philosophically random 
>effects as random no matter what, and leave them in the model even 
>if they don't appear to be statistically significant', and given 
>that with small numbers of random-effect levels this approach is 
>likely to lead to numerical difficulties in most (??) mixed model 
>packages (warnings, errors, or low estimates of the variance), what 
>should one do?
><snip>


================================================================
Robert A. LaBudde, PhD, PAS, Dpl. ACAFS  e-mail: ral at lcfltd.com
Least Cost Formulations, Ltd.            URL: http://lcfltd.com/
824 Timberlake Drive                     Tel: 757-467-0954
Virginia Beach, VA 23464-3239            Fax: 757-467-2947

"Vere scire est per causas scire"