# STAT425HW7Solution - STAT425HW7Solution Wenjing Yin...

• Homework Help
• 9
• 100% (10) 10 out of 10 people found this document helpful

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

The preview shows page 2 - 4 out of 9 pages.
STAT425HW7SolutionWenjing YinNovember 19, 2016setwd("~/Desktop/FALL 2016/STAT 425")library(faraway)## Warning: packagefarawaywas built under R version 3.2.3#### Attaching package:faraway## The following object is masked _by_.GlobalEnv :####choccakelibrary(lme4)## Warning: packagelme4was built under R version 3.2.5## Loading required package: Matrix## Warning: packageMatrixwas built under R version 3.2.5Problem 1 (14 pts)(a)video =read.table("videoviews.dat",header =TRUE)mod =lm(log(Views)~Channel,data =video)anova(mod)## Analysis of Variance Table#### Response: log(Views)##Df Sum Sq Mean Sq F valuePr(>F)## Channel5 30.1516.03023.241 0.02241 *## Residuals 24 44.6551.8606## ---## Signif. codes:0***0.001**0.01*0.05.0.11If Channel were a fixed factor, the p-value for Channel effect is 0.02241.(b)video.reml <-lmer(log(Views) ~ (1|Channel),data=video)summary(video.reml)1
## Linear mixed model fit by REML [ lmerMod ]## Formula: log(Views) ~ (1 | Channel)##Data: video#### REML criterion at convergence: 109.6#### Scaled residuals:##Min1QMedian3QMax## -1.3711 -0.7462 -0.14520.55742.6275#### Random effects:##GroupsNameVariance Std.Dev.##Channel(Intercept) 0.83390.9132##Residual1.86061.3640## Number of obs: 30, groups:Channel, 6#### Fixed effects:##Estimate Std. Error t value## (Intercept)4.05970.44839.055(c)sigma.alpha.hat.sq =0.9132^2sigma.hat.sq =1.3640^2lambda = sigma.alpha.hat.sq/(sigma.alpha.hat.sq+sigma.hat.sq)lambda## [1] 0.309503(d)video.ml <-lmer(log(Views) ~ (1|Channel),data=video,REML=FALSE)sigma.alpha.hat.sq =0.7956^2sigma.alpha.hat.sq## [1] 0.6329794(e)nullmod <-lm(log(Views) ~1,data=video)llrts <-as.numeric(2* (logLik(video.ml) -logLik(nullmod)))llrts# the log-likelihood-ratio test statistic## [1] 2.822395pchisq(llrts,1,lower=FALSE)# ordinary LRT p-value## [1] 0.09295759(f)2
set.seed(1)lrstats <-numeric(10000)for(i in1:10000){y <-unlist(simulate(nullmod))nullsim <-lm(y ~1)altsim <-lmer(y ~ (1|Channel),data=video,REML=FALSE)lrstats[i] <-as.numeric(2* (logLik(altsim) -logLik(nullsim)))}pval <-mean(lrstats >= llrts)pval# parametrically bootstrapped LRT p-value## [1] 0.022se.pval <-sqrt(pval*(1-pval)/10000)se.pval## [1] 0.001466833(g)ranef(video.reml)## \$Channel##(Intercept)

Course Hero member to access this document

Course Hero member to access this document

End of preview. Want to read all 9 pages?