Assignment 10 — Solution

Question 1

Read in and subset the data

# read in data
> srrs2 <- read.table("http://www.unc.edu/courses/2007spring/enst/562/001/assignments/assign10/radon.txt", header=T, sep=",")
# select only Minnesota data
> minn<-srrs2[srrs2$state=="MN",]
# obtain 6-character string to identify counties
> minn$countycode<-substr(as.character(minn$county),1,6)
# log-transform radon while adjusting zero values if necessary
> minn$logradon <- log(ifelse(minn$activity==0, .05, minn$activity))

• Create the graph.
  • As noted in the hints we need to calculate means for each value of floor and then connect the means by a line segment. This is done inside the panel function with tapply(y,x,mean). It's important to use the generic variables y and x inside the panel function. This causes the means to be obtained separately for each county instead of for the entire data set.
  • The function panel.lines then connects the two means (if there are two means) using the x-coordinates 1 and 2 in each panel.

# set background to white for older versions of R
> trellis.par.set(col.whitebg())
> xyplot(logradon~factor(floor)|countycode,data=minn, layout=c(17,5), strip=strip.custom(par.strip.text=list(cex=0.6)), panel=function(x,y) {
panel.xyplot(x,y)
mymeans<-tapply(y,x,mean)
panel.lines(as.numeric(names(mymeans))+1,mymeans,col=2)
}, xlab='Floor (0=basement, 1=first floor)', ylab='log(radon)')

Fig. 1 Lattice graph of household log radon levels by county for lowest floor in house

Question 2

• The unconditional means model is a model with no predictors, just an intercept at each level. In conventional multilevel model notation the unconditional means model is the following.

where

The composite equation for this model is obtained by inserting the level-2 equation into the level-1 equation.

• The lme function of R uses the composite version of the model in its syntax. The fixed and random parts are specified separately.

> library(nlme)
> lme(fixed=logradon~1,random=~1|county,data=minn,method='ML')->out0
> out0
Linear mixed-effects model fit by maximum likelihood
Data: minn
Log-likelihood: -1139.845
Fixed: logradon ~ 1
(Intercept)
1.308058

Random effects:
Formula: ~1 | county
      (Intercept)  Residual
StdDev: 0.3056031 0.8090579

Number of Observations: 919
Number of Groups: 85

• The useful information from this model is contained in the variance components. We can extract those with the VarCorr function.

> VarCorr(out0)
county = pdLogChol(1)
              Variance    StdDev
(Intercept) 0.09339322 0.3056031
Residual    0.65457471 0.8090579

From the output we learn that σ² = 0.655 and τ² = 0.093. Thus there is far more variability of household log radon levels within counties than between counties. This implies that the structure in this data set may not be particularly important. It also tells us that in terms of modeling we should invest most of our effort in trying to "explain" variability within counties rather than between counties, because that's where most of the variability lies. Thus level-1 predictors should prove to be more valuable than level-2 predictors.

• The intraclass correlation coefficient is

and can be calculated in R as follows.

> as.numeric(VarCorr(out0)[1,1])/ (as.numeric(VarCorr(out0)[1,1])+ as.numeric(VarCorr(out0)[2,1]))
[1] 0.1248626

• The correlation is not very large again suggesting that the households within a county are not particularly similar. This is probably not too surprising given that the county "structure" was not part of the original sampling design, but just a level added to the mix after the data were collected.

Question 3

Part 1

• Adding floor as a level-1 predictor to the unconditional means model yields the following multilevel model.

with the same probability models for the random effects as before. The corresponding composite equation is the following.

• The composite equation provides the format for fitting the model in R.

> lme(fixed=logradon~1+floor, random=~1|county, data=minn, method='ML')->out1
> summary(out1)
Linear mixed-effects model fit by maximum likelihood
Data: minn
     AIC      BIC    logLik
2195.602 2214.896 -1093.801

Random effects:
Formula: ~1 | county
       (Intercept) Residual
StdDev:  0.3230887 0.765721

Fixed effects: logradon ~ 1 + floor
                Value Std.Error  DF   t-value p-value
(Intercept) 1.4589422 0.0515057 833 28.325837       0
floor      -0.7036136 0.0713706 833 -9.858592       0
Correlation:
      (Intr)
floor -0.292

Standardized Within-Group Residuals:
        Min          Q1        Med         Q3        Max
-4.88589873 -0.60768113 0.01106789 0.62978548 3.44052086

Number of Observations: 919
Number of Groups: 85

• One piece of evidence suggesting that floor should be retained in the model is the statistical test provided by the row labeled "floor" in the summary table. The t-test for the hypothesis,

is highly significant (the reported p-value is rounded to 0) leading us to conclude that there is a statistically significant linear relationship between log radon and floor. In truth because floor is a dichotomous variable coded 0 and 1, we can also conclude that mean log radon levels are different in houses with basements compared to houses without basements.

• A second piece of evidence that this model should be preferred is obtained by comparing the AIC of this model to that of the unconditional means model.

> sapply(list(out0,out1),AIC)
[1] 2285.690 2195.602

• The AIC of the current model is much lower than that of the unconditional means model. Since these models only differ in the presence of the floor variable we can conclude that floor is an important predictor of log radon levels.

Part 2

• To assess whether intercepts should be random or fixed, we need to compare this model to a model with a fixed intercept, the model we've called the common pooling model. In this model we ignore the structure of the data (counties) and just fit a single regression model with floor as the predictor.

> out1a<-lm(logradon~floor,data=minn)

• There is an extension of the ordinary likelihood ratio test that can be used to compare these two models (since they are nested models), but we didn't cover it. The reason an extension is needed is because the hypothesized value of zero for the variance component lies on the boundary of the parameter space. We can use AIC to compare them though.

> sapply(list(out1,out1a),AIC)
[1] 2195.602 2273.466

The AIC is dramatically lower in the current model suggesting that allowing the intercepts to be random is a major improvement.

Part 3

• Because the predictor that we added, floor, is a level-1 predictor the appropriate R² to use is one which looks at the change in the estimated level-1 variance, , from a model with and a model without the predictor floor. This is a comparison of the current model with the unconditional means model.

> VarCorr(out0)
county = pdLogChol(1)
              Variance    StdDev
(Intercept) 0.09339322 0.3056031
Residual    0.65457471 0.8090579
> VarCorr(out1)
county = pdLogChol(1)
             Variance    StdDev
(Intercept) 0.1043863 0.3230887
Residual    0.5863287 0.7657210

• The appropriate comparison is between the entries labeled "Residual" in the two sets of output. Formally,

• We can calculate this quantity in R by copying and pasting the values on the screen or more formally as follows.

> (as.numeric(VarCorr(out0)[2,1])- as.numeric(VarCorr(out1)[2,1]))/ as.numeric(VarCorr(out0)[2,1])
[1] 0.1042601

• So the addition of the floor variable has contributed to a 10% reduction in between household variability.

Question 4

• The random slopes and intercepts model is the following multilevel model.

The probability models for the random effects change because we've added an additional random effect.

The corresponding composite equation is the following.

• The composite equation provides the format for fitting the model in R.

> lme(fixed=logradon~1+floor, random=~1+floor|county, data=minn, method='ML')->out2

• Models with multiple random effects can be difficult to fit in the frequentist paradigm so it's always a good idea to verify that model actually converged to a solution. The variance components can be informative in this regard.

> VarCorr(out2)
county = pdLogChol(1 + floor)
             Variance    StdDev  Corr
(Intercept) 0.1162887 0.3410114 (Intr)
floor       0.1120414 0.3347258 -0.319
Residual    0.5726597 0.7567428

• When there are convergence problems one of the variance components is often returned with a value that is approximately zero. That's not the case here.
• There is an extension of the likelihood ratio test for comparing these two models, but we didn't study it. Instead I compare these modes using AIC.

> sapply(list(out1,out2),AIC)
[1] 2195.602 2197.084

As we can see the AIC of the random slopes and intercepts model actually went up relative to random intercepts model. Therefore we should prefer the simpler model and stick with fixed slopes but random intercepts.

Question 5

• I load the county data, subset it to include just the Minnesota counties, and merge it with the household data. For the merge I use the by.x and by.y arguments because the key fields, cntyfips and ctfips, have different names in the two files and there are other fields with the same name in both files.

> cty <- read.table ("http://www.unc.edu/courses/2007spring/enst/562/001/assignments/assign10/cty.txt", header=T, sep=",")
> minn2<-cty[cty$st=='MN',]
> merge(minn,minn2,by.x="cntyfips",by.y="ctfips")->merge3

• The variable log(Uppm) varies at the county level and therefore is added as a predictor to the level-2 equations. Technically we have a choice, we can add it to the level-2 equation for the intercept or the level-2 equation for the slope. I was not very clear in my instructions on this point (my bad). But "adding it to the random intercepts model" turns out to mean "add it to the level-2 equation for the intercept" as we'll now see.
• Here's the level-2 equations with log(Uppm) as a predictor of the intercept.

where as before . The composite equation is the following.

• We fit this in R as follows.

> lme(fixed=logradon~1+floor+log(Uppm), random=~1|county, data=merge3, method='ML')->out3
> summary(out3)
Linear mixed-effects model fit by maximum likelihood
Data: merge3
     AIC      BIC    logLik
2157.913 2182.029 -1073.956

Random effects:
Formula: ~1 | county
       (Intercept) Residual
StdDev: 0.1457421  0.7689218

Fixed effects: logradon ~ 1 + floor + log(Uppm)
                Value  Std.Error  DF  t-value p-value
(Intercept) 1.4627406 0.03756813 833 38.93568       0
floor      -0.6788426 0.06969627 833 -9.74001       0
log(Uppm)   0.7180305 0.09064038  83  7.92175       0
Correlation:
         (Intr)  floor
floor    -0.362
log(Uppm) 0.156 -0.011

Standardized Within-Group Residuals:
        Min          Q1        Med         Q3        Max
-5.12860182 -0.60172541 0.03181403 0.64952352 3.32903088

Number of Observations: 919
Number of Groups: 85

• The row of the summary table labeled "log(Uppm)" is a test of

As we can see from the reported p-value, this test is highly significant. Thus we should retain this variable. We can also reach this conclusion by comparing the AIC of this model to that of a model without log(Uppm).

> sapply(list(out1,out3),AIC)
[1] 2195.602 2157.913

The addition of log(Uppm) has led to a dramatic decrease in AIC. Thus we should prefer a model with log(Uppm) as a predictor.

• Since log(Uppm) is a level-2 predictor, the appropriate pseudo-R² to examine is one that looks at the change in level-2 variance.

Here the baseline model is the random intercepts model without the log(Uppm) predictor. I extract the components I need from the VarCorr object and carry out the calculations.

> VarCorr(out1)
county = pdLogChol(1)
             Variance    StdDev
(Intercept) 0.1043863 0.3230887
Residual    0.5863287 0.7657210
> VarCorr(out3)
county = pdLogChol(1)
              Variance    StdDev
(Intercept) 0.02124076 0.1457421
Residual    0.59124069 0.7689218
> (as.numeric(VarCorr(out1)[1,1])- as.numeric(VarCorr(out3)[1,1]))/ as.numeric(VarCorr(out1)[1,1])
[1] 0.7965177

• The calculated R² is gigantic. 80% of the variability in the model intercepts between counties is explained by their linear relationship with average county soil uranium levels

Course Home Page