[R-sig-ME] Random vs. fixed effects
Robert A LaBudde
ral at lcfltd.com
Sat Apr 24 07:17:19 CEST 2010
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"
More information about the R-sig-mixed-models
mailing list