srrs2 <- read.table ("srrs2.dat", header=T, sep=",")
mn <- srrs2$state=="MN"
radon <- srrs2$activity[mn]
log.radon <- log (ifelse (radon==0, .1, radon))
floor <- srrs2$floor[mn] # 0 for basement, 1 for first floor
n <- length(radon)
y <- log.radon
x <- floor
county.name <- as.vector(srrs2$county[mn])
uniq <- unique(county.name)
J <- length(uniq)
county <- rep (NA, J)
for (i in 1:J){
county[county.name==uniq[i]] <- i
}###pooled model
lm.pooled<-lm(formula=y~x)
summary(lm.pooled)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6293 -0.5383 0.0342 0.5603 2.5486
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.32674 0.02972 44.640 <2e-16 ***
## x -0.61339 0.07284 -8.421 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8226 on 917 degrees of freedom
## Multiple R-squared: 0.07178, Adjusted R-squared: 0.07077
## F-statistic: 70.91 on 1 and 917 DF, p-value: < 2.2e-16###unpooled model
lm.unpooled<-lm(formula=y~x+factor(county)-1)
summary(lm.unpooled)##
## Call:
## lm(formula = y ~ x + factor(county) - 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.14595 -0.45405 0.00065 0.45376 2.65987
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## x -0.72054 0.07352 -9.800 < 2e-16 ***
## factor(county)1 0.84054 0.37866 2.220 0.026701 *
## factor(county)2 0.87482 0.10498 8.333 3.23e-16 ***
## factor(county)3 1.52870 0.43946 3.479 0.000530 ***
## factor(county)4 1.55272 0.28897 5.373 1.00e-07 ***
## factor(county)5 1.43257 0.37866 3.783 0.000166 ***
## factor(county)6 1.51301 0.43672 3.464 0.000558 ***
## factor(county)7 2.01216 0.20243 9.940 < 2e-16 ***
## factor(county)8 1.98958 0.37999 5.236 2.08e-07 ***
## factor(county)9 1.00304 0.23931 4.191 3.07e-05 ***
## factor(county)10 1.56391 0.31099 5.029 6.04e-07 ***
## factor(county)11 1.40113 0.33828 4.142 3.80e-05 ***
## factor(county)12 1.73025 0.37821 4.575 5.49e-06 ***
## factor(county)13 1.03872 0.30881 3.364 0.000804 ***
## factor(county)14 1.98838 0.20325 9.783 < 2e-16 ***
## factor(county)15 1.33797 0.37999 3.521 0.000453 ***
## factor(county)16 0.66486 0.53487 1.243 0.214204
## factor(county)17 1.27480 0.38221 3.335 0.000890 ***
## factor(county)18 1.12155 0.21913 5.118 3.83e-07 ***
## factor(county)19 1.33831 0.09541 14.026 < 2e-16 ***
## factor(county)20 1.80032 0.43672 4.122 4.13e-05 ***
## factor(county)21 1.73399 0.25227 6.873 1.23e-11 ***
## factor(county)22 0.63679 0.30905 2.060 0.039663 *
## factor(county)23 1.39999 0.53613 2.611 0.009183 **
## factor(county)24 2.10162 0.25267 8.318 3.64e-16 ***
## factor(county)25 1.95072 0.20243 9.636 < 2e-16 ***
## factor(county)26 1.36058 0.07422 18.332 < 2e-16 ***
## factor(county)27 1.77336 0.30978 5.725 1.45e-08 ***
## factor(county)28 1.24159 0.34115 3.639 0.000290 ***
## factor(county)29 1.05600 0.43672 2.418 0.015818 *
## factor(county)30 0.92576 0.22807 4.059 5.39e-05 ***
## factor(county)31 2.02057 0.33828 5.973 3.45e-09 ***
## factor(county)32 1.23629 0.37821 3.269 0.001124 **
## factor(county)33 2.06187 0.37821 5.452 6.58e-08 ***
## factor(county)34 1.59044 0.43946 3.619 0.000314 ***
## factor(county)35 0.81920 0.28897 2.835 0.004695 **
## factor(county)36 2.95897 0.53613 5.519 4.55e-08 ***
## factor(county)37 0.40209 0.25227 1.594 0.111345
## factor(county)38 1.86772 0.37999 4.915 1.07e-06 ***
## factor(county)39 1.74807 0.33860 5.163 3.05e-07 ***
## factor(county)40 2.31580 0.37866 6.116 1.48e-09 ***
## factor(county)41 1.96715 0.26759 7.351 4.69e-13 ***
## factor(county)42 1.36098 0.75642 1.799 0.072343 .
## factor(county)43 1.60224 0.25543 6.273 5.69e-10 ***
## factor(county)44 1.04099 0.28609 3.639 0.000291 ***
## factor(county)45 1.29541 0.21101 6.139 1.28e-09 ***
## factor(county)46 1.21461 0.33828 3.591 0.000349 ***
## factor(county)47 0.88393 0.53613 1.649 0.099583 .
## factor(county)48 1.14812 0.25227 4.551 6.13e-06 ***
## factor(county)49 1.70211 0.21010 8.102 1.93e-15 ***
## factor(county)50 2.49321 0.75642 3.296 0.001022 **
## factor(county)51 2.16504 0.37821 5.724 1.45e-08 ***
## factor(county)52 1.92769 0.43672 4.414 1.15e-05 ***
## factor(county)53 1.25080 0.43741 2.860 0.004348 **
## factor(county)54 1.30676 0.15802 8.270 5.28e-16 ***
## factor(county)55 1.61799 0.26885 6.018 2.64e-09 ***
## factor(county)56 1.10110 0.43946 2.506 0.012415 *
## factor(county)57 0.76218 0.30905 2.466 0.013855 *
## factor(county)58 1.86092 0.37866 4.915 1.07e-06 ***
## factor(county)59 1.72178 0.37999 4.531 6.73e-06 ***
## factor(county)60 1.27939 0.53487 2.392 0.016979 *
## factor(county)61 1.15873 0.13389 8.654 < 2e-16 ***
## factor(county)62 1.98301 0.33860 5.856 6.80e-09 ***
## factor(county)63 1.67070 0.43741 3.820 0.000144 ***
## factor(county)64 1.84784 0.22817 8.099 1.97e-15 ***
## factor(county)65 1.29912 0.53487 2.429 0.015357 *
## factor(county)66 1.66574 0.20648 8.067 2.50e-15 ***
## factor(county)67 1.80312 0.21101 8.545 < 2e-16 ***
## factor(county)68 1.09002 0.26743 4.076 5.02e-05 ***
## factor(county)69 1.24245 0.37821 3.285 0.001062 **
## factor(county)70 0.86763 0.07096 12.227 < 2e-16 ***
## factor(county)71 1.49184 0.15174 9.832 < 2e-16 ***
## factor(county)72 1.57990 0.23920 6.605 7.08e-11 ***
## factor(county)73 1.79176 0.53487 3.350 0.000845 ***
## factor(county)74 0.98704 0.37821 2.610 0.009223 **
## factor(county)75 1.72372 0.43741 3.941 8.80e-05 ***
## factor(county)76 2.00844 0.37866 5.304 1.45e-07 ***
## factor(county)77 1.82168 0.28609 6.367 3.17e-10 ***
## factor(county)78 1.28569 0.33956 3.786 0.000164 ***
## factor(county)79 0.61488 0.37866 1.624 0.104785
## factor(county)80 1.32952 0.11181 11.890 < 2e-16 ***
## factor(county)81 2.70953 0.43946 6.166 1.09e-09 ***
## factor(county)82 2.23001 0.75642 2.948 0.003286 **
## factor(county)83 1.62292 0.21048 7.711 3.57e-14 ***
## factor(county)84 1.64535 0.20987 7.840 1.38e-14 ***
## factor(county)85 1.18652 0.53487 2.218 0.026801 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7564 on 833 degrees of freedom
## Multiple R-squared: 0.7671, Adjusted R-squared: 0.7431
## F-statistic: 31.91 on 86 and 833 DF, p-value: < 2.2e-16library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack#basic MLM with just random intercept for county
M1<-lmer(y~x+(1|county),REML=FALSE)
summary(M1)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: y ~ x + (1 | county)
##
## AIC BIC logLik deviance df.resid
## 2171.7 2190.9 -1081.8 2163.7 915
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4071 -0.6164 0.0056 0.6398 3.4288
##
## Random effects:
## Groups Name Variance Std.Dev.
## county (Intercept) 0.1053 0.3245
## Residual 0.5703 0.7552
## Number of obs: 919, groups: county, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.46116 0.05124 28.516
## x -0.69264 0.07036 -9.844
##
## Correlation of Fixed Effects:
## (Intr)
## x -0.290Here we pull in county-level data on uranium concentrations in the soil.
srrs2.fips <- srrs2$stfips*1000 + srrs2$cntyfips
cty <- read.table ("cty.dat", header=T, sep=",")
usa.fips <- 1000*cty[,"stfips"] + cty[,"ctfips"]
usa.rows <- match (unique(srrs2.fips[mn]), usa.fips)
uranium <- cty[usa.rows,"Uppm"]
u <- log (uranium)
u.full<-u[county]# Extended model adding uranium but still only letting county effect be random
M2<-lmer(y~x+u.full+(1|county),REML=FALSE)
summary(M2)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: y ~ x + u.full + (1 | county)
##
## AIC BIC logLik deviance df.resid
## 2132.8 2156.9 -1061.4 2122.8 914
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9976 -0.6163 0.0307 0.6561 3.3794
##
## Random effects:
## Groups Name Variance Std.Dev.
## county (Intercept) 0.02127 0.1458
## Residual 0.57499 0.7583
## Number of obs: 919, groups: county, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.46427 0.03714 39.421
## x -0.66644 0.06865 -9.708
## u.full 0.72320 0.08965 8.067
##
## Correlation of Fixed Effects:
## (Intr) x
## x -0.361
## u.full 0.154 -0.011You can add the fixed effects and random effects intercept estimates to get the county-specific estimates.
#county-specific intercepts (already summed)
coef(M2)## $county
## (Intercept) x u.full
## 1 1.446360 -0.6664446 0.7232005
## 2 1.477562 -0.6664446 0.7232005
## 3 1.475294 -0.6664446 0.7232005
## 4 1.564351 -0.6664446 0.7232005
## 5 1.471756 -0.6664446 0.7232005
## 6 1.441176 -0.6664446 0.7232005
## 7 1.581533 -0.6664446 0.7232005
## 8 1.502609 -0.6664446 0.7232005
## 9 1.403162 -0.6664446 0.7232005
## 10 1.464858 -0.6664446 0.7232005
## 11 1.523062 -0.6664446 0.7232005
## 12 1.473036 -0.6664446 0.7232005
## 13 1.483560 -0.6664446 0.7232005
## 14 1.552980 -0.6664446 0.7232005
## 15 1.450082 -0.6664446 0.7232005
## 16 1.434359 -0.6664446 0.7232005
## 17 1.402969 -0.6664446 0.7232005
## 18 1.495567 -0.6664446 0.7232005
## 19 1.385946 -0.6664446 0.7232005
## 20 1.478770 -0.6664446 0.7232005
## 21 1.502003 -0.6664446 0.7232005
## 22 1.273595 -0.6664446 0.7232005
## 23 1.437307 -0.6664446 0.7232005
## 24 1.579944 -0.6664446 0.7232005
## 25 1.579079 -0.6664446 0.7232005
## 26 1.432816 -0.6664446 0.7232005
## 27 1.450991 -0.6664446 0.7232005
## 28 1.469662 -0.6664446 0.7232005
## 29 1.477798 -0.6664446 0.7232005
## 30 1.447268 -0.6664446 0.7232005
## 31 1.516209 -0.6664446 0.7232005
## 32 1.439860 -0.6664446 0.7232005
## 33 1.531057 -0.6664446 0.7232005
## 34 1.473829 -0.6664446 0.7232005
## 35 1.456401 -0.6664446 0.7232005
## 36 1.549861 -0.6664446 0.7232005
## 37 1.322404 -0.6664446 0.7232005
## 38 1.576315 -0.6664446 0.7232005
## 39 1.484922 -0.6664446 0.7232005
## 40 1.530801 -0.6664446 0.7232005
## 41 1.512364 -0.6664446 0.7232005
## 42 1.456722 -0.6664446 0.7232005
## 43 1.488732 -0.6664446 0.7232005
## 44 1.350955 -0.6664446 0.7232005
## 45 1.371060 -0.6664446 0.7232005
## 46 1.422598 -0.6664446 0.7232005
## 47 1.432893 -0.6664446 0.7232005
## 48 1.400647 -0.6664446 0.7232005
## 49 1.477518 -0.6664446 0.7232005
## 50 1.490677 -0.6664446 0.7232005
## 51 1.531461 -0.6664446 0.7232005
## 52 1.481265 -0.6664446 0.7232005
## 53 1.421987 -0.6664446 0.7232005
## 54 1.307791 -0.6664446 0.7232005
## 55 1.528542 -0.6664446 0.7232005
## 56 1.429738 -0.6664446 0.7232005
## 57 1.356566 -0.6664446 0.7232005
## 58 1.469035 -0.6664446 0.7232005
## 59 1.469165 -0.6664446 0.7232005
## 60 1.437526 -0.6664446 0.7232005
## 61 1.459781 -0.6664446 0.7232005
## 62 1.502191 -0.6664446 0.7232005
## 63 1.455592 -0.6664446 0.7232005
## 64 1.533374 -0.6664446 0.7232005
## 65 1.426588 -0.6664446 0.7232005
## 66 1.574772 -0.6664446 0.7232005
## 67 1.554066 -0.6664446 0.7232005
## 68 1.491669 -0.6664446 0.7232005
## 69 1.413569 -0.6664446 0.7232005
## 70 1.252753 -0.6664446 0.7232005
## 71 1.432915 -0.6664446 0.7232005
## 72 1.442803 -0.6664446 0.7232005
## 73 1.463372 -0.6664446 0.7232005
## 74 1.373305 -0.6664446 0.7232005
## 75 1.491761 -0.6664446 0.7232005
## 76 1.486276 -0.6664446 0.7232005
## 77 1.513860 -0.6664446 0.7232005
## 78 1.508870 -0.6664446 0.7232005
## 79 1.333256 -0.6664446 0.7232005
## 80 1.442874 -0.6664446 0.7232005
## 81 1.571802 -0.6664446 0.7232005
## 82 1.485489 -0.6664446 0.7232005
## 83 1.402955 -0.6664446 0.7232005
## 84 1.542854 -0.6664446 0.7232005
## 85 1.427437 -0.6664446 0.7232005
##
## attr(,"class")
## [1] "coef.mer"#fixed effects
fixef(M2)## (Intercept) x u.full
## 1.4642651 -0.6664446 0.7232005#random effects
ranef(M2)## $county
## (Intercept)
## 1 -0.0179049553
## 2 0.0132967914
## 3 0.0110287558
## 4 0.1000862363
## 5 0.0074905000
## 6 -0.0230891610
## 7 0.1172682332
## 8 0.0383440827
## 9 -0.0611033042
## 10 0.0005928621
## 11 0.0587970219
## 12 0.0087706144
## 13 0.0192948986
## 14 0.0887146888
## 15 -0.0141831962
## 16 -0.0299061999
## 17 -0.0612961327
## 18 0.0313021006
## 19 -0.0783189236
## 20 0.0145052057
## 21 0.0377374237
## 22 -0.1906697471
## 23 -0.0269580812
## 24 0.1156793235
## 25 0.1148136002
## 26 -0.0314492806
## 27 -0.0132743951
## 28 0.0053968655
## 29 0.0135324582
## 30 -0.0169968356
## 31 0.0519434864
## 32 -0.0244053538
## 33 0.0667917505
## 34 0.0095636906
## 35 -0.0078641968
## 36 0.0855960132
## 37 -0.1418611800
## 38 0.1120497108
## 39 0.0206570459
## 40 0.0665354713
## 41 0.0480987057
## 42 -0.0075429376
## 43 0.0244664122
## 44 -0.1133097791
## 45 -0.0932052603
## 46 -0.0416670319
## 47 -0.0313718174
## 48 -0.0636185948
## 49 0.0132526872
## 50 0.0264114516
## 51 0.0671959196
## 52 0.0169994943
## 53 -0.0422783415
## 54 -0.1564744846
## 55 0.0642764988
## 56 -0.0345269336
## 57 -0.1076987708
## 58 0.0047696890
## 59 0.0048996570
## 60 -0.0267387480
## 61 -0.0044843852
## 62 0.0379256899
## 63 -0.0086726566
## 64 0.0691085347
## 65 -0.0376774205
## 66 0.1105065369
## 67 0.0898010353
## 68 0.0274040943
## 69 -0.0506964759
## 70 -0.2115118076
## 71 -0.0313500327
## 72 -0.0214618586
## 73 -0.0008935427
## 74 -0.0909600970
## 75 0.0274960635
## 76 0.0220112930
## 77 0.0495946579
## 78 0.0446046828
## 79 -0.1310090098
## 80 -0.0213916011
## 81 0.1075369683
## 82 0.0212238364
## 83 -0.0613104471
## 84 0.0785885116
## 85 -0.0368282736
##
## with conditional variances for "county"We can also allow variation across county in the association between floor of measurement and radon concentration.
#changing optimizer from default due to convergence issues
M3<-lmer(y~x+u.full+(1+x|county),REML=FALSE,control = lmerControl(optimizer ="Nelder_Mead"))
summary(M3)## Linear mixed model fit by maximum likelihood ['lmerMod']
## Formula: y ~ x + u.full + (1 + x | county)
## Control: lmerControl(optimizer = "Nelder_Mead")
##
## AIC BIC logLik deviance df.resid
## 2131.6 2165.4 -1058.8 2117.6 912
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1425 -0.6099 0.0321 0.6374 3.4975
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## county (Intercept) 0.01398 0.1182
## x 0.11909 0.3451 0.26
## Residual 0.56032 0.7485
## Number of obs: 919, groups: county, 85
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.46115 0.03480 41.982
## x -0.64069 0.08620 -7.432
## u.full 0.77057 0.08679 8.879
##
## Correlation of Fixed Effects:
## (Intr) x
## x -0.257
## u.full 0.195 -0.050