STAT425HW7Solution - STAT425HW7Solution Wenjing Yin...

STAT425HW7Solution
Wenjing Yin
November 19, 2016
setwd ( "~/Desktop/FALL 2016/STAT 425" )
library (faraway)
## Warning: package faraway was built under R version 3.2.3
##
## Attaching package: faraway
## The following object is masked _by_ .GlobalEnv :
##
## choccake
library (lme4)
## Warning: package lme4 was built under R version 3.2.5
## Loading required package: Matrix
## Warning: package Matrix was built under R version 3.2.5
Problem 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 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 1
If 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:
## 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.055
(c)
sigma.alpha.hat.sq = 0.9132 ^ 2
sigma.hat.sq = 1.3640 ^ 2
lambda = 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 ^ 2
sigma.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.822395
pchisq (llrts, 1 , lower= FALSE ) # ordinary LRT p-value
## [1] 0.09295759
(f)
2

set.seed ( 1 )
lrstats <- numeric ( 10000 )
for(i in 1 : 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.022
se.pval <- sqrt (pval*( 1 -pval)/ 10000 )
se.pval
## [1] 0.001466833
(g)
ranef (video.reml)
## $Channel
## (Intercept)