########## # Problem 11.1 ########## It is OK if student consider other models. Just make sure their answer are consistent. I tried many models and only provide answer for the model with weeks, treat, and their interaction effect in the fixed-effect components. The random components only contain treat effect. (b) Read my output. Note that I used zero sum constraint. My fitted model is wt=54.20857-1.32857*I(treat=control)+3.45143*I(treat=thiouracil)-2.12286*I(treat=thyroxine) +(23.57762-2.90238*I(treat=control)-6.46762*I(treat=thiouracil)+3.56524*I(treat=thyroxine))*weeks +gamma_{0,subject}+gamma_{1,subject}*weeks+N(0, 4.347983^2), where (gamma_{0,subject},gamma_{1,subject}) follows bivariate normal with mean zero, variances 5.700354^2 and 3.760148^2, correlation -0.133. Keep in mind to interpret the i, ii, iii. The fixed-effect is the mean of the change. The randon-effect is the variation of the change. (c) I use the likelihood ratio test. The test statistic is 33.47 based on 10 degrees of freedom. Therefore, it is significant. You can also use the bootstrap method. The bootstrap p-value is 0. Since the test is for fixed-effect, it is not necsseay to use the bootstrap method. Note: You may use the bootstrap method to test significance of random effect, but not required here. It is OK if you use lmer function. ######### # Problem 11.2 ######### (b) Only ajwtr1 is not significant. The coefficient for time means that the housping price decreasing in time, which is not reasonable. (d) ajwtr1 is still insignificant. The coefficient for time is positive, which is reasonable. (e) MSA are related to metropolitan statistical areas. It can be used as a random effect. I omit the interpretation. Make sure to interpret the fixed effect intercept as the mean of intercept and the random effect as its variation. The rest components does not have variations in the model. (g) I got -5.965818 loglikelihood ratio statistic value. Clearly, it is not significant. You may also use the bootstrap method. Clearly it is not needed. ######### # R output ######### > library(grid) > library(lattice) > library(nlme) > library(faraway) > ######## > > ratdring <- read.table("u:\\stat526\\dataset\\ascdata\\ratdrink.txt",h=T) > summary(ratdring) wt weeks subject treat Min. : 46.0 Min. :0 Min. : 1 control :50 1st Qu.: 71.0 1st Qu.:1 1st Qu.: 7 thiouracil:50 Median :100.0 Median :2 Median :14 thyroxine :35 Mean :100.8 Mean :2 Mean :14 3rd Qu.:122.5 3rd Qu.:3 3rd Qu.:21 Max. :189.0 Max. :4 Max. :27 > > data(ratdring) Warning message: In data(ratdring) : data set ‘ratdring’ not found > summary(ratdring) wt weeks subject treat Min. : 46.0 Min. :0 Min. : 1 control :50 1st Qu.: 71.0 1st Qu.:1 1st Qu.: 7 thiouracil:50 Median :100.0 Median :2 Median :14 thyroxine :35 Mean :100.8 Mean :2 Mean :14 3rd Qu.:122.5 3rd Qu.:3 3rd Qu.:21 Max. :189.0 Max. :4 Max. :27 > > # (b) > > options(contrasts=c("contr.sum", "contr.sum")) > > mod.b.1 <- lme(wt~weeks*treat,random=~weeks|subject,data=ratdring) > summary(mod.b.1) Linear mixed-effects model fit by REML Data: ratdring AIC BIC logLik 903.0698 931.668 -441.5349 Random effects: Formula: ~weeks | subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 5.700354 (Intr) weeks 3.760148 -0.133 Residual 4.347983 Fixed effects: wt ~ weeks * treat Value Std.Error DF t-value p-value (Intercept) 54.20857 1.2922776 105 41.94808 0.0000 weeks 23.57762 0.7814324 105 30.17231 0.0000 treat1 -1.32857 1.7695240 24 -0.75081 0.4601 treat2 3.45143 1.7695240 24 1.95048 0.0629 weeks:treat1 2.90238 1.0700204 105 2.71245 0.0078 weeks:treat2 -6.46762 1.0700204 105 -6.04439 0.0000 Correlation: (Intr) weeks treat1 treat2 wks:t1 weeks -0.250 treat1 -0.091 0.023 treat2 -0.091 0.023 -0.400 weeks:treat1 0.023 -0.091 -0.250 0.100 weeks:treat2 0.023 -0.091 0.100 -0.250 -0.400 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.83135768 -0.54990678 0.04002798 0.58230401 2.03659601 Number of Observations: 135 Number of Groups: 27 > > # (c) > > gg <- lm(wt~weeks,data=ratdring) > xx <- model.matrix(gg) > > > mod.b.2 <- lme(wt~weeks,random=~weeks|subject,data=ratdring) > summary(mod.b.2) Linear mixed-effects model fit by REML Data: ratdring AIC BIC logLik 928.5403 945.8824 -458.2702 Random effects: Formula: ~weeks | subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 5.961342 (Intr) weeks 5.957620 -0.404 Residual 4.347989 Fixed effects: wt ~ weeks Value Std.Error DF t-value p-value (Intercept) 54.44444 1.317694 107 41.31796 0 weeks 23.18148 1.176683 107 19.70070 0 Correlation: (Intr) weeks -0.433 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.83279170 -0.55110985 -0.01505971 0.55847226 2.05813684 Number of Observations: 135 Number of Groups: 27 > Like.obs <- 2*(logLik(mod.b.1)-logLik(mod.b.2)) > Like.obs 'log Lik.' 33.47047 (df=10) > > BETA <- mod.b.2$coefficients$fixed > > sigmasq.mu <- 5.961342^2 > sigmasq.weeks <- 5.957620^2 > rho.mu.weeks <- -0.404 > sigma <- 4.347989 > > > AA <- cbind(c(sigmasq.mu,rho.mu.weeks*sqrt(sigmasq.mu*sigmasq.weeks)),c(rho.mu.weeks*sqrt(sigmasq.mu*sigmasq.weeks),3.760148^2)) > AA [,1] [,2] [1,] 35.53760 -14.34823 [2,] -14.34823 14.13871 > > > library(MASS) > BB <- mvrnorm(27,mu=rep(0,2),Sigma=AA) > random.mu <- rep(BB[,1],rep(5,27)) > random.weeks <- rep(BB[,2],rep(5,27)) > yy <- xx%*%BETA+random.mu+random.weeks*ratdring$weeks+rnorm(135,mean=0,sd=sigma) > > > lrstat <- rep(0,1000) > for(i in 1:1000){ + AA <- cbind(c(sigmasq.mu,rho.mu.weeks*sqrt(sigmasq.mu*sigmasq.weeks)),c(rho.mu.weeks*sqrt(sigmasq.mu*sigmasq.weeks),3.760148^2)) + BB <- mvrnorm(27,mu=rep(0,2),Sigma=AA) + random.mu <- rep(BB[,1],rep(5,27)) + random.weeks <- rep(BB[,2],rep(5,27)) + yy <- xx%*%BETA+random.mu+random.weeks*ratdring$weeks+rnorm(135,mean=0,sd=sigma) + mod.a <- lme(yy~weeks*treat,random=~weeks|subject,data=ratdring) + mod.n <- lme(yy~weeks,random=~weeks|subject,data=ratdring) + lrstat[i] <- 2*(logLik(mod.a)-logLik(mod.n)) + } > > pp <- mean(lrstat>Like.obs) > pp [1] 0 > sqrt(pp*(1-pp)/1000) [1] 0 ########### # Problem 2 ########### > data(hprice) > summary(hprice) narsp ypc perypc regtest rcdum ajwtr msa time Min. :3.920 Min. :12535 Min. :-2.054 Min. :13.00 0:279 0:189 1 : 9 Min. :1 1st Qu.:4.264 1st Qu.:16609 1st Qu.: 3.535 1st Qu.:18.00 1: 45 1:135 2 : 9 1st Qu.:3 Median :4.412 Median :18454 Median : 3.964 Median :20.00 3 : 9 Median :5 Mean :4.484 Mean :18769 Mean : 4.268 Mean :20.42 4 : 9 Mean :5 3rd Qu.:4.575 3rd Qu.:20323 3rd Qu.: 5.711 3rd Qu.:22.00 5 : 9 3rd Qu.:7 Max. :5.563 Max. :33383 Max. : 8.788 Max. :29.00 6 : 9 Max. :9 (Other):270 > > # (b) > > mod.1 <- lm(narsp~ypc+perypc+regtest+rcdum+ajwtr+time,data=hprice) > summary(mod.1) Call: lm(formula = narsp ~ ypc + perypc + regtest + rcdum + ajwtr + time, data = hprice) Residuals: Min 1Q Median 3Q Max -0.31386 -0.10810 -0.01525 0.08547 0.55594 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.765e+00 9.398e-02 29.425 < 2e-16 *** ypc 7.029e-05 4.358e-06 16.128 < 2e-16 *** perypc -1.372e-02 5.074e-03 -2.704 0.007216 ** regtest 2.954e-02 3.103e-03 9.520 < 2e-16 *** rcdum1 -7.440e-02 1.618e-02 -4.599 6.15e-06 *** ajwtr1 -1.797e-02 1.001e-02 -1.796 0.073482 . time -1.767e-02 5.128e-03 -3.445 0.000647 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1651 on 317 degrees of freedom Multiple R-squared: 0.7572, Adjusted R-squared: 0.7526 F-statistic: 164.7 on 6 and 317 DF, p-value: < 2.2e-16 > # (d) > > hprice$ypc.first <- rep(hprice$ypc[9*(0:35)+1],rep(9,36)) > > mod.2 <- lm(narsp~ypc.first+perypc+regtest+rcdum+ajwtr+time,data=hprice) > summary(mod.2) Call: lm(formula = narsp ~ ypc.first + perypc + regtest + rcdum + ajwtr + time, data = hprice) Residuals: Min 1Q Median 3Q Max -0.37278 -0.10546 -0.01540 0.07898 0.53906 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.395e+00 1.060e-01 22.606 < 2e-16 *** ypc.first 8.865e-05 5.275e-06 16.806 < 2e-16 *** perypc -8.128e-03 4.971e-03 -1.635 0.103 regtest 3.014e-02 3.036e-03 9.929 < 2e-16 *** rcdum1 -7.560e-02 1.581e-02 -4.782 2.66e-06 *** ajwtr1 -1.927e-02 9.795e-03 -1.967 0.050 . time 3.711e-02 3.785e-03 9.803 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1619 on 317 degrees of freedom Multiple R-squared: 0.7662, Adjusted R-squared: 0.7618 F-statistic: 173.2 on 6 and 317 DF, p-value: < 2.2e-16 > # (e) > > mod.3 <- lme(narsp~ypc.first+perypc+regtest+rcdum+ajwtr+time,random=~1|msa,data=hprice) > summary(mod.3) Linear mixed-effects model fit by REML Data: hprice AIC BIC logLik -560.4889 -526.6588 289.2444 Random effects: Formula: ~1 | msa (Intercept) Residual StdDev: 0.1536575 0.0738323 Fixed effects: narsp ~ ypc.first + perypc + regtest + rcdum + ajwtr + time Value Std.Error DF t-value p-value (Intercept) 2.4009031 0.29064017 286 8.260741 0.0000 ypc.first 0.0000886 0.00001521 31 5.829574 0.0000 perypc -0.0091479 0.00229776 286 -3.981214 0.0001 regtest 0.0301616 0.00874728 31 3.448109 0.0016 rcdum1 -0.0756850 0.04555680 31 -1.661332 0.1067 ajwtr1 -0.0192643 0.02823530 31 -0.682275 0.5001 time 0.0368038 0.00172932 286 21.282175 0.0000 Correlation: (Intr) ypc.fr perypc regtst rcdum1 ajwtr1 ypc.first -0.748 perypc -0.045 0.002 regtest -0.508 -0.174 -0.005 rcdum1 -0.571 0.332 0.004 0.310 ajwtr1 -0.169 0.182 -0.001 0.040 -0.178 time -0.043 0.001 0.395 -0.002 0.002 0.000 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.38462421 -0.59377571 -0.02501864 0.57411299 2.85776453 Number of Observations: 324 Number of Groups: 36 > # (g) > > mod.4 <- lme(narsp~ypc.first+perypc+regtest+time,random=~1|msa,data=hprice) > summary(mod.4) Linear mixed-effects model fit by REML Data: hprice AIC BIC logLik -570.4547 -544.0984 292.2273 Random effects: Formula: ~1 | msa (Intercept) Residual StdDev: 0.1577741 0.07383233 Fixed effects: narsp ~ ypc.first + perypc + regtest + time Value Std.Error DF t-value p-value (Intercept) 2.0458563 0.23068646 286 8.868558 0e+00 ypc.first 0.0001007 0.00001422 33 7.084191 0e+00 perypc -0.0091335 0.00229778 286 -3.974946 1e-04 regtest 0.0355055 0.00848977 33 4.182154 2e-04 time 0.0368080 0.00172933 286 21.284596 0e+00 Correlation: (Intr) ypc.fr perypc regtst ypc.first -0.698 perypc -0.053 0.000 regtest -0.416 -0.349 -0.006 time -0.053 0.000 0.395 -0.002 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.37925756 -0.59579797 -0.01311158 0.58008158 2.84995440 Number of Observations: 324 Number of Groups: 36 > > Like.obs <- 2*(logLik(mod.3)-logLik(mod.4)) > Like.obs 'log Lik.' -5.965818 (df=9)