# STAT425_HW_7_Yilun_Luo.docx - STAT425 HW#7 Yilun Luo...

• Homework Help
• 9
• 50% (2) 1 out of 2 people found this document helpful

This preview shows page 1 - 4 out of 9 pages.

STAT425 HW#7Yilun LuoDecember 2, 20161part(a)videoviews =read.table("videoviews.dat",header =TRUE)fit1 =lm(log(Views) ~Channel, data =videoviews)anova(fit1)## Analysis of Variance Table## ## Response: log(Views)## Df Sum Sq Mean Sq F value Pr(>F) ## Channel 5 30.151 6.0302 3.241 0.02241 *## Residuals 24 44.655 1.8606 ## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The p-value for channel effect is 0.02241.part(b)library(lme4)## Loading required package: Matrixfit2 =lmer(log(Views) ~(1|Channel), data =videoviews)summary(fit2)## Linear mixed model fit by REML ['lmerMod']## Formula: log(Views) ~ (1 | Channel)## Data: videoviews## ## REML criterion at convergence: 109.6## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -1.3711 -0.7462 -0.1452 0.5574 2.6275 ## ## Random effects:## Groups Name Variance Std.Dev.## Channel (Intercept) 0.8339 0.9132 ## Residual 1.8606 1.3640 ## Number of obs: 30, groups: Channel, 6##
## Fixed effects:## Estimate Std. Error t value## (Intercept) 4.0597 0.4483 9.055part(c)Intraclasscorrelation =0.8339/(0.8339+1.8606)Intraclasscorrelation## [1] 0.3094823Intra-class correlation = 0.3094823.part(d)fit3 =lmer(log(Views) ~(1|Channel), data =videoviews, REML =FALSE)summary(fit3)## Linear mixed model fit by maximum likelihood ['lmerMod']## Formula: log(Views) ~ (1 | Channel)## Data: videoviews## ## AIC BIC logLik deviance df.resid ## 115.7 119.9 -54.9 109.7 27 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -1.4236 -0.7363 -0.1520 0.5881 2.5750 ## ## Random effects:## Groups Name Variance Std.Dev.## Channel (Intercept) 0.6329 0.7956 ## Residual 1.8606 1.3640 ## Number of obs: 30, groups: Channel, 6## ## Fixed effects:## Estimate Std. Error t value## (Intercept) 4.0597 0.4093 9.919part(e)fit4 =lm(log(Views) ~1, data =videoviews)llrts =as.numeric(2*(logLik(fit3) -logLik(fit4)))llrts## [1] 2.822395pchisq(llrts, 1, lower =FALSE)## [1] 0.09295759part(f)llrts1 =numeric(10000)for(i in 1:10000){ y =unlist(simulate(fit4))
nullsim =lm(y ~1,data =videoviews)altsim =lmer(y ~(1|Channel), data =videoviews, REML =FALSE)llrts1[i] =as.numeric(2*(logLik(altsim) -logLik( nullsim))) }## Warning in optwrap(optimizer, devfun, getStart(start, rho\$lower, rho\$pp),:## convergence code 3 from bobyqa: bobyqa -- a trust region step failed to## reduce qp.val =mean(llrts1 >=llrts)p.val## [1] 0.0219