Saturday, October 27, 2007

Mixed Models - random coefficient models

library(nlme)

coplot(Y ~ EP|No, data=petrol)

Petrol <- petrol
Petrol[,2:5] <- scale(Petrol[,2:5],scale=F)

# Y~EP for each No; no overal mean
pet1.lm <- lm(Y~No/EP-1, data=Petrol)

# Y~ mu_i + beta*EP; no overal mean
pet2.lm <- lm(Y~No-1+EP, data=Petrol)
pet3.lm <- lm(Y~SG+VP+V10+EP, data=Petrol)

# Random intercept model
pet3.lme <- lme(Y~SG+VP+V10+EP, random=~1|No, data=Petrol)

# To compare likelihood, use "ML"
pet3.lme <- update(pet3.lme, method="ML")
anova(pet3.lme, pet3.lm)
pet4.lme <- update(pet3.lme, fixed=Y~V10+EP)
# pet4.lme <- lme(Y~V10+EP, random=~1|No, method="ML", data=Petrol)
anova(pet3.lme, pet4.lme)

## Random intercept + Random coefficient
pet5.lme <- lme(Y~V10+EP, random=~1+EP|No, method="ML", data=Petrol)
anova(pet4.lme, pet5.lme)

No comments: