1 #!/usr/local/bin/R
  2 
  3 # Xavier Fernández i Marín
  4 # December 2004
  5 # http://xavier-fim.net/
  6 
  7 # reply the analysis done with HLM
  8 
  9 # load the example datasets from HLM/AppendixA/Examples
 10 # Download them directly. Requires Internet connection
 11 hsb1 <- read.table(file="http://xavier-fim.net/R/scripts/hsb1.dat")
 12 hsb2 <- read.table(file="http://xavier-fim.net/R/scripts/hsb2.dat")
 13 
 14 # get the variables and basic information of each dataset
 15 names(hsb1)
 16 length(hsb1$ID)
 17 summary(as.data.frame(hsb1))
 18 
 19 names(hsb2)
 20 length(hsb2$ID)
 21 summary(as.data.frame(hsb2))
 22 
 23 # load the package to do Linear Mixed-Effects Models
 24 library(nlme)
 25 
 26 # create a new dataset with the merged data
 27 ds <- merge(as.data.frame(hsb1), as.data.frame(hsb2))
 28 
 29 ### model specification
 30 
 31 # 1
 32 mod.1 <- lme(fixed = MATHACH ~ 1,
 33              random = ~ 1 | ID,
 34              data=ds, method="ML"
 35              )
 36 
 37 # 2
 38 mod.2 <- lme(fixed = MATHACH ~ SES,
 39              random = ~ 1 | ID,
 40              data=ds, method="ML"
 41              )
 42 
 43 # 3
 44 mod.3 <- lme(fixed = MATHACH ~ SES + SECTOR,
 45              random = ~ 1 | ID,
 46              data=ds, method="ML"
 47              )
 48 
 49 # 4
 50 mod.4 <- lme(fixed = MATHACH ~ SES + SECTOR,
 51              random = ~ 1 + SES | ID,
 52              data=ds, method="ML"
 53              )
 54 
 55 # 5
 56 mod.5 <- lme(fixed = MATHACH ~ SES * SECTOR,
 57              random = ~ 1 + SES | ID,
 58              data=ds, method="ML"
 59              )
 60 
 61 
 62 # Extract information about the models
 63 print(mod.1)      # basic information
 64 summary(mod.1)    # complete information
 65 summary(mod.2)
 66 summary(mod.3)
 67 summary(mod.4)
 68 summary(mod.5)
 69 
 70 # attributes of an 'lmeObject'
 71 mod.1$dims     # dimensions 
 72 mod.1$coef     # coefficients
 73 mod.1$call     # the formula used to produce the object
 74 mod.1$residuals     # residuals of the model
 75 mod.1$varFix   # matrix with the fixed effects
 76 mod.1$sigma    # estimated within-group error standard deviation
 77 
 78 anova(mod.1)
 79 
 80 
 81 # Plots
 82 plot(mod.1)
 83      # by default, the plot of an lme object is 
 84      # form=resid(., type = "p") ~ fitted(.)
 85      # that corresponds to the standardized residuals vs predicted values
 86 
 87 plot(mod.1, resid(.) ~ fitted(.) | SECTOR)
 88      # idem, but comparing the sectors 
 89 
 90 qqnorm(mod.1)     # normality of the residuals
 91 
 92 
 93 
 94 
 95 interaction.plot(ds$SES, ds$SECTOR, ds$MATHACH, col=c("blue", "red"))
 96      # show one of the plotting abilities of R
 97 
 98 # standardized residuals vs predicted values, by sector
 99 plot(mod.3, resid(., type = "p") ~ fitted(.) | SECTOR, abline = 0)
100 
101 # boxplots with the residuals by school 
102 #  (it's just an example, there are too many objects in level 2)
103 plot(mod.4, ID ~ resid(.))
104 
105