Homework 3 Solutions - Multilevel Hierarchical Models HW03 Solution Problem One library(lme4 library(arm mn.radon < read.table\"mn-radon attach(mn.radon

Homework 3 Solutions - Multilevel Hierarchical Models HW03...

This preview shows page 1 out of 10 pages.

You've reached the end of your free preview.

Want to read all 10 pages?

Unformatted text preview: Multilevel & Hierarchical Models HW03 Solution November 17, 2013 Problem One library(lme4) library(arm) mn.radon <- read.table("mn-radon.txt") attach(mn.radon) n <- length(radon) j <- length(unique(county)) y <- log.radon x <- floor u.full <- log.uranium M0 <- lmer( y ~ 1 + (1|county)) M1 <- lmer(y ~ x + (1|county)) M2 <- lmer(y ~ x + u.full + (1|county)) M3 <- lmer(y ~ x + u.full +(1+ x|county)) anova(M0,M1,M2,M3) Data: Models: M0: y ~ 1 + (1 | county) M1: y ~ x + (1 | county) M2: y ~ x + u.full + (1 | county) M3: y ~ x + u.full + (1 + x | county) Df AIC BIC logLik Chisq Chi Df Pr(>Chisq) M0 3 2261.2 2275.7 -1127.6 M1 4 2171.7 2190.9 -1081.8 91.5876 1 < 2.2e-16 *** M2 5 2132.9 2157.0 -1061.4 40.7947 1 1.691e-10 *** M3 7 2131.7 2165.4 -1058.8 5.2002 2 0.07426 . display(M0) number of obs: 919, AIC = 2265.4, DIC = deviance = 2255.2 display(M1) number of obs: 919, AIC = 2179.3, DIC = deviance = 2163.7 display(M2) groups: county, 85 2251 groups: county, 85 2156 1 number of obs: 919, AIC = 2144.2, DIC = deviance = 2122.9 display(M3) number of obs: 919, AIC = 2142.6, DIC = deviance = 2117.7 groups: county, 85 2111.5 groups: county, 85 2106.7 code for producing all the residuals plot for M2, M3 is similar: par(mfrow = c(2,3)) plot(yhat.marg(M2),r.marg(M2),xlab="marginal fitted values",ylab="marginal residuals") abline(h=0,lwd = 2, col = "red3") plot(yhat.cond(M2),r.cond(M2),xlab="conditional fitted values",ylab="conditional residuals") abline(h=0,lwd = 2, col = "red3") plot(yhat.reff(M2),r.reff(M2),xlab="random effect fitted values",ylab="random effect residuals") abline(h=0,lwd = 2, col = "red3") n2<-qqnorm(r.marg(M2), main="qqplot for Marginal Residuals") abline(lm(y~x,data=n2),col="red",lwd=2) n2<-qqnorm(r.cond(M2), main="qqplot for Conditional Residuals") abline(lm(y~x,data=n2),col="red",lwd=2) n2<-qqnorm(r.reff(M2), main="qqplot for Marginal Residuals") abline(lm(y~x,data=n2),col="red",lwd=2) xyplot(r.reff(M2)~yhat.reff(M2)|as.factor(county)) We summarize the results in the following table. Table 1: Model Selection Criteria for Radon Models AIC BIC DIC M0 2261 2276 2251 M1 2172 2191 2156 M2 2133 2157 2112 M3 2132 2165 2107 We can see that based on AIC, M2 and M3 are equally the best models. based on BIC, M2 is the best model. based on DIC, M3 is the best model. Plot for M2: 2 1.0 qq 0.5 q 1.0 0.00 q q q q q qqqq qq q q qq q q q q qqq q q q q qq q q q qq q q q qq q q q q q qq q q qqqq qqq q qq q q q 1.5 qq qqq qq q q q qq q q q qqq q q q qq q q qqqqqq q q q qqq q q q qq q q q q qq qq qq q q qq q q qq qqq qqqq q qq q q q qqq q q q q qq qq q q qqqqqqqqqq qq q q qqqqq q q qq q q q q q q q qq q q q q q qq q qq q q qq q qq q q q q q qq q q q q qq q qq q qq qq qq q q q q q q qqq qq qq qq qqq q qq qq q qq q qq q q q qq q q q qq qq qqqqq qq q q q q q q q qqqqqqqqq qqqq q qqqqq q qqqqqq q qqqq qqqqq q q q q qqq qq q qq qq q q q q qq q q qq q q qq qqqqq qq qq qq q q qq q q qqq qqq qq qqq q q qq q qq qq q qq qqq qqqqqqq qq qq q q q q qq qqq qqqqq q q q q q q qqq q qq q q q q qq q q qqq q q qqq q qq qqqqqqqqq qq qqqqqqqq q q qqq qqqqq qq q q qq q q qqqqqqqqqqq qq q q q q qq qq q q q qq q q qq q q q q qqqqqq q q q q qqq qq q q qq q q qq q q q q qqq qqqq qqqq q q q qqq q q q qq qqq q q q qq q qq q qq 0.10 q −0.10 q random effect residuals 2 −1 0 1 q 0.5 −2 −3 q −4 q q q q qq qq q q q q q qq q q q qq qq q q q q q q qq qq q q q qq q qq q q q qq q q q q q qq qq qq q q q qq q q q qq q q q q qq q qq qq q q q q qq qq q qq qq q q q q qq q q qq qqqq q qq q q q qq q q qqq qq q q qq qq q qq q q qqq qqqq q q q q q qqq q qq qqqq q q q q q q q q q q q q q qq qq q qq q q qq q q q qqqq q q q q qq qq q q qq q q q q qq qq q q q q q qqq q q qqq q q qq q q q qq q qq q q q qq qq q q q qq q q q q q q q qq q q q q qqqq q q q qq q q qq qq q q q q q qq q q qq q q q qqq qq qq qqqq qqq q qq qq q q qqq q qq qq q q qqq q q q q qq q qq q q q q q qq q q q qqqqq q q qq qqq qqq q qq q q q qq q q q q q q q q q q qq q qq qq q qq q qqqqq q qq q qq qq q qq q q q q q qqqq qq qq qq qqqq q q q q qq q qq q qqq q q q q qqq q q qq q q q qq q q q q qqq q q q q q q q q qq qq qq q qq q q q q qq q q q q q q q q qq q q qq q qq qq q qq q q q q q q q q q q qq q q q q q qq q q qq q q q qq q qqq qq q q q qq qqq q q q q qq q q q qq q q q q q qq q q qq q q q qq q q qqq q q qq q q q qq q q q qq q qq qq qq q q q q q q q qq q q q q −0.20 q q q −4 −2 q q conditional residuals 1 2 q q q q q −1 0 q −3 marginal residuals q q q q qq q qq q qq qq q q qq q q q q q qqq q qq q q q q qq q q q qq q q q q q qq q q q qq q qqq q q qq q q q qqq qq q q q q qq q q q qq q q q qq q q q qqq q qq q qq q q qq qq q qq qq q q qqq q q q q qqq q qq qq q qq q q qq q qq qq q q q q q qq q q qq q q qqqqq qq qq qqq qq q q q q q q qqqq q q q q q q q qq q qq q q qq qq q q qq qq q qq qqqq q q q q qq q q q q qq q q q q qqq q q q qq q q qq q q q q qq q q q q qqqq qq q q q q q qqq q q q q q qq q q q qqqq q qq q q q q q qq qqq q q q q qq q q q qq q q q q qq q q qq q q q q q q q qqqqq q q qqqq qqq qq q qq q qq q qq qq q qq qq q q q qq qqq qq qq qqq qq q q q q qq q q q qqq q q q q qq q q q qq q q qq qq q q qq q q qq q q qqqqq q q q q qqqqq q q q q q qq qqqq q q q qqq qqq q q q q q q q q q qq q qq q q qq q q qq q q q q q qq q q q q qq q q q q q q qq q q qq q q q q qqq qq q q qqq qqq q q q q q q q q q qq qq q q q qqq q q q q q q qq qq q q q qq qq q q qq q q q qq q q q q q q q q qqq qqq q q q q qq q qq q q q q q qq q q q q q q q q qq q q qq q q q q q q q q q q q qq 1.5 q q qqqqqq qq qqq q qqqq qqqqq qq q q q qqqqqqqq q q q qq q q q qq q −2 −1 0 1 2 q 3 4 marginal fitted values conditional fitted values random effect fitted values qqplot for Marginal Residuals qqplot for Conditional Residuals qqplot for Marginal Residuals −1 0 1 2 Theoretical Quantiles 3 0.10 0.00 −0.10 Sample Quantiles q q q q −3 −1 0 1 2 Theoretical Quantiles 3 q qqq q qqq qq q q qqq q q q qq qqq qq qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq qq q q q q q q q q q q q q q q q q q qq qq qq qq qq qq q q q q q q q qq q q q q q q q q q q q q q q q q −0.20 2 1 0 −1 −2 −4 q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −3 q q q −3 q q Sample Quantiles 0 −1 −2 −3 −4 Sample Quantiles 1 2 q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q 3 q q q qqqqqqqq qqqqqqqq qqqqqqq qqqqqqq q qqqq qqq −3 −1 0 1 2 Theoretical Quantiles 3 −2 0 2 4 81 0.1 0.0 −0.1 −0.2 82 qq q 83 84 85 qqq qqq q q q q q q qqq q qq qq 71 72 qqq qq qq q qqqq qqqq qqq qq 61 0.1 0.0 −0.1 −0.2 −2 0 2 4 62 q qq qq qqqq qqqq qqq qq q 73 74 75 76 qq q q q 77 q qq q q q qq q q qq q q 63 52 qq q q 80 qqqq q qqqq qqqq qq qq qq q q 64 65 66 67 qqq qq q q q qq q qq qq q qq q q qq qq qq q qq 68 69 q q q 41 42 qq q q q q 53 54 55 56 57 qq q q q q q qq 43 70 qq qq q 32 33 qq q qq qq q 44 qqq q 34 58 qqq q q q 59 60 qq qq 0.1 0.0 −0.1 −0.2 q qq qq q 45 46 q q q qq q q q q qq q 31 q qq qq qq qq q 35 47 q 36 q 25 50 q qq qq q q 38 qq q q 39 40 qqq q qq q q qq q q q 24 49 q qq q q qq qq q 37 qq q q q q qq q 48 qqqq q q 21 0.1 0.0 −0.1 −0.2 22 qq q q q q 23 qqqq qq q qqq q qq qq q q q q 11 qq qq 1 0.1 0.0 −0.1 −0.2 0.1 0.0 −0.1 −0.2 qq q q qqq qqqq qqq qq r.reff(M2) 79 qqqqqq q qqqqq qqqqq qqq qqq 51 0.1 0.0 −0.1 −0.2 78 qq q 26 −2 0 2 4 qq qq qqqqq q qqqqq qq qq qqqq qq q 28 29 qq qq q q q q 18 30 19 qq q qq qq qq qq q 12 13 qq q qqq qq q 2 3 14 15 16 q qqqq qqqq qqq qqq qq qq qq −2 0 2 4 17 qqqq qqq q q qq q 4 q q 5 6 q qq q qq q 27 0.1 0.0 −0.1 −0.2 qq q qq q 7 qqq qq q q yhat.reff(M2) Plot for M3: 4 q q qqqq q qqqqq qqq qq q qq q 8 9 qq qq qq q −2 0 2 4 20 q qq qq qq q qqq qq qq −2 0 2 4 10 q −2 0 2 4 q qq q q 0.1 0.0 −0.1 −0.2 q q 0.2 0.0 −0.2 q q q q q q q q qq qq q qq q qq q qq qq q q q qq qq qq qq q q q qq qq q qq qqq q q qq qq q q q q q q qq qq qqq q q q q q q qq q qqq qq q q q qqq q qq q qq qq qq qqq q qq qqqq qq qq q qq qqq qq q qq q q qqqq q qq qq qqq q q q q q q q qq q q q q q q qq q qq q q qq qqqq q qq q qq q qq qqq qq qq qqqqq q q qq qq q q q q q qq qq qq q q q q q qqqq qqq q q q qqq q q qq q q q qq qq q q qq qq q qq qqqqqq q qq qqqqq qqq q q q qq q qq q q q q q q q qq q q q qqq q q q qqq qq qq q q qq q q qq q q q qqqq qq q q q qq q q q q qq q qq q qq q qq qq q q q q q qq qq q q q qq qq q q q qq q qqq q q q qq q q qqqq q q q q qqqqqqq q qq qqq q q q qq qq qq q q q qq q qq q q qqq q q qqqq q qq q qq qq q q qq q qq qq q qqq q q q qq q q q q q qqq q q q qq q qqq q q qq q q qq qq q q q qq qq q qq q qqqq q q qqqqqq q q q q qq qq q q q qq qq qq q q qq q q q q qq q qq q q q qq q q q qq qq q q qq q q q q qq q qq q q qq q q qq q q q q q qq q q q qqq q q qq q q q q q q q q random effect residuals 2 1 q q qq q qq q q q q qq q q q q q q q q q q q −3 q q q q q q q q 0 qq qq qq qqq q q q qq q q q qqq q qq q q q qq q q q qqqqq q qq q qq q qq q q q qq q q qq q q q qq q q q q qq q q q q qqq qqq q q qq q qq q q q qq q q qqq qqq qq q qq q q q qq q q q qq qq q q qq q qqq qq q q q qqqq qqq q q q qq q q q qq qq q qq qq q q qq qqq q q q qq q q qq q q q qq q q q qq q q q q qq q q q q qq q qq q q q qqqq q q q qq q q qq q q qq q q qqq qq q q q qq q q q qq q qq qq q q q q qq q qq qq q qq q qq q q q qq q −1 q qq q q qq qq q q q qq q q qq q q q qq q q q q qq qq q q q qq q q qqqq q qqqq qq qq q q qq q q qqqqq qq qq qqq q qq qq q q q q qq q q qq q qq q qqq q qq q qq qqqq q qqq qqq q q q qq q q q qqq q q qq q q q q q qq qq q q qq qqq q qq qqqqq q q qqq q q qq qqq q q q qq qqq q qq q q q q qqqq qqq q qq qq q qqq qqq q qq qqq qq qq q qqqq q q q qq q qqq qq q q q qqq qqq qq qqq q q qqq q q q q qq qqq q q q q qq q q qqqqq qq qq q qq q q qq qq q q q qq q q q qq q q q q qqq q qq qq q q q q q qq q q qqq qqq q qq qq q q qqq qq q qq q q 0.4 q q q −0.4 q q −2 qq q q qq q qq qq qq q q q q q qq q q q q conditional residuals 2 1 0 −1 −2 q −3 marginal residuals q qq q q q q q qq q qq q q qq q q qq q q qq q qq q q qq qq q q q q qqq q q q q q qq q q q q qq qq q q qq q qq q qq qqq q q q q qq q q qq qqq qqq q q qq q qqq qqq qqq q qqq qq q q qq q qqqqqqqqqq qq qqq qq qqq q q q q qqqqqqqqqqqqq q q q q qqqqqqqqqqqq q q q q qqqqqqqqq q q qqqqq qqqq q q qqqqqq q q qq qq qqq qqqqqqqqqq q q q qqqq q qq qqqqqqqqqqqqqqq qq qqqqqqqq qq q qqqqqqqq qqqq qqqqqqqqqqqqqqqq q qqqqqqqq qq q q q qqqqqqqqqqqq qq q qqqqqqqq qq q q q qq qq q q q qqqqqqqqqqqqq qqqqq qqq q qq q qqqqqq q q qq q q q qqqqq q qq q qq qqq q q qqqq qqqqqqqqq q q qqq qqq qqqq qqq qqq qqqq q q qq qqqq q qq q q q qq q q q qq q q q qq q q q q qq qq qq q q q q q q qq q qq qqqq q qq qq q q q qq q q qq q q qqqq q q qq qq q q qq q qq qq q qqq qqqqqqqqq q q q q q q q qqqqqqq q q q q q qq qqqqq qq q q q qqqqq q q qq q q q qqq q q q q q q qq q qqq q q q q q q 0.5 1.0 q −4 1.5 −0.6 −4 q q 0.0 0.5 1.0 1.5 qq q q qq qq q qq qq qq −2 −1 0 q 1 2 3 4 marginal fitted values conditional fitted values random effect fitted values qqplot for Marginal Residuals qqplot for Conditional Residuals qqplot for Marginal Residuals 0.4 qq q q qq qq q q q q q −0.2 0.0 0.2 q q q q qq q q q q q q q q q q q qq q q qq qq qq qq qq qq q qq qq qq q qq qq qq qq qq qq q q q q q qq qqq qqq qqq qqq qqq qqq q q q q q q q q q qq q q q q q q q q q q qqq qqq qqq qqq qqq qq q q q q qq qq q q q q q q q qq q q qq qq −0.4 q q q Sample Quantiles 2 1 0 −1 −2 q q q q q −3 Sample Quantiles 0 −1 −2 −3 qq q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −3 −4 q −1 0 1 2 Theoretical Quantiles 3 −0.6 q q q −4 Sample Quantiles 1 2 q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q −3 −1 0 1 2 Theoretical Quantiles 5 3 qq q −3 −1 0 1 2 Theoretical Quantiles 3 −2 01234 81 0.2 −0.2 −0.6 82 q q q q 71 72 qqq qq qq q 61 62 q qqq q qq q qqqq qqq qqq qq q 51 52 r.reff(M3) qq q q q q q 41 0.2 −0.2 −0.6 42 q qq q q q 31 32 qq qq q qqq q 21 0.2 −0.2 −0.6 83 qq qq qq 84 85 qqq qqq q q q q qq q qqqq qqq qqq q q qq q q 0.2 −0.2 −0.6 −2 01234 qq qq q q 22 q 11 qq qq qq q q q 12 qq q 73 q q 63 qq q 53 qq q 43 q q qq q q q 33 qq q 23 q q 13 qq q qq q 74 75 qq qq 64 68 qq qq q qq q 58 59 qq q qq q q q q q 34 35 q qq 25 qqq q q qq q 14 qqq q q qq qq 47 48 q qq qq q q 37 38 q qqq q q 26 27 qq qq qqqqq q qqqq qqqq qqqq qq qqqqq qq qq q 16 qq q 17 q q 1 qq q q −2 01234 2 qqqq qqqq qqqq qqq qq q qq q 3 qq −2 01234 4 q qq q q 5 6 q qq qq q −2 01234 q qq 7 qq qqq qq q q −2 01234 q q qq 50 q 40 q qq q q 28 q q q 18 qqq q 29 qqq q q 19 qq q qqq qq q 8 qq qq 30 0.2 −0.2 −0.6 qq q qq qq 20 q q qqqqq qqqqq qqq qq q qq q qq qq 9 q qqq qq q 0.2 −0.2 −0.6 q 39 q q q 15 60 qqq q qq q q q q qq q q 0.2 −0.2 −0.6 qqqqqq q qqqqq qqqqq qqqq qqq 49 36 q q qqq q qq q q 70 qq qq qq q q qq q 0.2 −0.2 −0.6 qqq q q q qq qq q q q 24 57 46 q q qqqq q qqqq qqqq qq qq qq qq 69 q qq q q qq qq q 56 45 80 qq q q 67 55 44 qq q 79 66 qq q qq q q qqqq qqq qq qq qq q qq q q q q 78 qq qq q qq q q qq 54 77 q qq q 65 qq qq q qq q q q qq q q 76 q qq 10 q 0.2 −0.2 −0.6 q qq q q −2 01234 yhat.reff(M3) All residual plots show similar properties between model2 and model3. Another way to do goodness of fit test is to construct certain test statistics and tell from the p-value. For example, if we want to compare M2 and M3, we can compare their deviance and degree of freedom and use chi-square test statistics to see if they are significantly different etc. We also produced some facet plots, but they were not very informative. We did not reach any clear interpretations for them, and didn’t find them useful. We show only the plots for M3 as examples. Below are facet plots by floor and by county (showing only the largest counties) for fits against residuals for the conditional, marginal, and random effects residuals. 6 Code used to produce plots is shown below. mod = M3 plotdata <- data.frame([email protected], fitted.re = fitted(M0), rmarg = r.marg(mod), rcond = r.cond(mod), rreff = r.reff(mod), ymarg = yhat.marg(mod), ycond = yhat.cond(mod), yreff = yhat.reff(mod)) plotdata = plotdata[plotdata$county %in% bigCounties,] qplot(data = plotdata, x = ymarg, y = rmarg, facets = ~ county) 7 qplot(data = plotdata, x = ycond, y = rcond, facets = ~ county) qplot(data = plotdata, x = yreff, y = rreff, facets = ~ county) Problem Two A L(λ) = ¯ e−nλ λnx x1 !x2 ! . . . xn ! The log likelihood is ¯ (λ) = log L(λ) = −nλ + nX log λ − log x1 !...xn ! ∂ (λ) ¯1 = −n + n X ∂λ λ set the first derivative as 0 and solve for λ, we have ˆ ¯ λ=X ˆ To verify that λ achieves the maximum, take second derivative of the log likelihood ∂ 2 (λ) ¯1 = −n X 2 < 0 ∂λ2 λ which means the log likelihood function is concave with respective to λ. B ¯1 From the last problem, we know the second derivative of the log likelihood function is −nX λ2 . Then E (− n ∂ 2 (λ) )=− ∂λ2 λ Then I (λ)−1 = ˆ λ n ˆ ˆ ¯ where λ is the MLE and we know that λ = X . So Iλ−1 = I (λ)−1 = ¯ X n. Thus the standard error is ¯ X n Problem Three We have likelihood and prior: L(p) ∝ px (1 − p)n−x and the prior distribution f (p) ∝ pα−1 (1 − p)β −1 Then the posterior distribution π (p) ∝ px+α−1 (1 − p)n−x+β −1 It is easy to see that the posterior distribution is also a Beta distribution with parameters α∗ = x + α, β ∗ = n − x + β . Given x = 5, n = 10, we try the following three values of α and β . 8 1. α = 2, β = 1 and α∗ = 7, β ∗ = 6. 2. α = 1, β = 2 and α∗ = 6, β ∗ = 7. 3. α = 3, β = 6 and α∗ = 8, β ∗ = 11. The density plots of these three distributions are in the following figure. x Different priors make the posterior different. For the first case the mean of prior is greater than n , while x it is less than n for the second and third cases. We can see the mean of posterior distribution is dominated by the prior distribution. 9 4 posterior distribution of p 2 1 0 pdf 3 alpha*=7,beta*=6 alpha*=6,beta*=7 alpha*=8,beta*=11 0.0 0.2 0.4 0.6 p 10 0.8 1.0 ...
View Full Document

  • Fall '13
  • BrianW.Junker

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture

  • Left Quote Icon

    Student Picture