# UPDATED 7 MAR 2006 TO INCLUDE INTERACTION IN MODEL 3 ############################################################################## # Fit the simulated examples in R # Read in data for Examples 1-3 y123 <- read.table ("y123.dat", header=TRUE) y <- y123\$y x <- y123\$x group <- y123\$group u <- read.table ("u.dat", header=TRUE)\$u ############################################################################## # Fit model 1 M1 <- lmer (y ~ x + (1 | group)) # Example 1: y_i ~ N (a_{group[i]} + b*x_i, sigma.y^2) # a_j ~ N (mu.a, sigma.a^2) # True (population) values from simulation: # mu.a=5, b=-3, sigma.a=2, sigma.y=6 # After fitting the model: estimates (+/- se's) from data: # mu.a=5.2(+/-1.8), b=-3.3(+/-0.7), sigma.a=3.2, sigma.y=6.3 ############################################################################## # Fit model 2 u.full <- u[group] M2 <- lmer (y ~ x + u.full + (1 | group)) # Example 2: y_i ~ N (a_{group[i]} + b*x_i, sigma.y^2) # a_j ~ N (gamma.0 + gamma.1*u, sigma.a^2) # True (population) values from simulation: # gamma.0=5, gamma.1=0, b=-3, sigma.a=2, sigma.y=6 # After fitting the model: estimates (+/- se's) from data: # gamma.0=1.2(+/-2.6), gamma.1=2.7(+/-1.3), b=-3.2(+/-0.7), sigma.a=2.6, sigma.y=6.3 ############################################################################## # Fit model 3 # UPDATED 7 MAR 2006 TO INCLUDE INTERACTION IN MODEL 3 M3 <- lmer (y ~ x + u.full + x:u.full + (1 + x | group)) # Example 3: y_i ~ N (a_{group[i]} + b_{group[i]}*x_i, sigma.y^2) # The J pairs (a_j,b_j) follow a bivariate normal distribution # with mean vector (gamma.0 + gamma.1*u[j], delta.0 + delta.1*u[j]) # and covariance matrix estimated from the data # True (population) values from simulation: # gamma.0=5, gamma.1=0, delta.0=-3, delta.1=0, sigma.a=2, sigma.b=0, rho.ab=0, sigma.y=6 # After fitting the model: estimates (+/- se's) from data: # gamma.0=-2.4(+/-5.2), gamma.1=5.0(+/-3.2), delta.0=-1.5(+/-2.0), delta.1=-1.1(+/-1.3), sigma.a=5.1, sigma.b=1.7, rho.ab=-0.9, sigma.y=6.1 ############################################################################## # Read in data for Example 4 y4 <- read.table ("y4.dat", header=TRUE) y <- y4\$y x <- y4\$x group <- y4\$group ############################################################################## # Fit model 4 M4 <- lmer (y ~ x + (1 | group), family=binomial(link="logit")) # Example 4: Pr(y_i=1) = invlogit (a_{group[i]} + b*x_i) # a_j ~ N (mu.a, sigma.a^2) # True (population) values from simulation: # mu.a=5, b=-3, sigma.a=2 # After fitting the model: estimates (+/- se's) from data: # mu.a=6.1(+/-1.5), b=-3.2(+/-0.7), sigma.a=2.5 ############################################################################## # Read in data for Example 5 y5 <- read.table ("y5.dat", header=TRUE) y <- y5\$y x <- y5\$x z <- y5\$z group <- y5\$group ############################################################################## # Fit model 5 log.z <- log(z) M5 <- lmer (y ~ x + (1 | group), offset=log.z, family=quasipoisson(link="log")) # Example 5: y_i ~ overdispersed Poisson (z_i*exp(a_{group[i]} + b*x_i)) # a_j ~ N (mu.a, sigma.a^2) # True (population) values from simulation: # mu.a=5, b=-3, sigma.a=2 # After fitting the model: estimates (+/- se's) from data: # mu.a=6.4(+/-0.5), b=-3.6(+/-0.1), sigma.a=1.7 ############################################################################## # Read in data for Example 6 y6 <- read.table ("y6.dat", header=TRUE) y <- y6\$y state <- y6\$state occupation <- y6\$occupation ############################################################################## # Fit model 6 state.occupation <- max(occupation)*(state - 1) + occupation M6 <- lmer (y ~ 1 + (1 | state) + (1 | occupation) + (1 | state.occupation)) # Example 6: y_i ~ N (mu + a_{state[i]} + b_{occupation[i]} + g_{state[i],occupation[i], sigma.y^2) # a_j ~ N (mu.a, sigma.a^2) # b_j ~ N (mu.b, sigma.b^2) # g_j ~ N (mu.g, sigma.g^2) # True (population) values from simulation: # mu.4, sigma.a=2, sigma.b=4, sigma.g=5, sigma.y=6 # After fitting the model: estimates (+/- se's) from data: # mu.4.5(+/-0.9), sigma.a=2.1, sigma.b=2.3, sigma.g=5.7, sigma.y=6.2