set.seed(100)
social.eff <- sample(ads.social,4)
ads.social <- 0:100/100
ads.normal <- 0:100/100
ads.animated <- 150:250/250
set.seed(100)
social.eff <- sample(ads.social,4)
normal.eff <- sample(ads.normal,4)
animated.eff <- sample(ads.animated,4)
##100 of each ad type
clicks <- c(runif(50*4)<social.eff, runif(50*4)<normal.eff, runif(50*4) <animated.eff) ##600 long; in
type <- as.factor(rep(c("A","B","C"),
each=200))             ## 600 long
ad <- as.factor(rep(rep(1:4,50),3))
ad2 <- as.factor(c(rep(1:4,50),rep(11:14,50), rep(21:24,50)))                           ##ad id; in blocks of type.
##
sub <- as.factor(rep(rep(1:50,each=4),3))
##Set the contrasts to be sum-to-zero so they are orthogonal
contrasts(type) <- contr.sum(3)
contrasts(ad) <- contr.sum(length(unique(ad)))
contrasts(ad2) <- contr.sum(length(unique(ad2)))
contrasts(sub) <- contr.sum(50)
ads <- data.frame(sub=sub,type=type,ad=ad,ad2=ad2,clicks=clicks)
ads[1:10,]
ads
table(ads$ad,ads$ad2)
pairwise.t.test(clicks,type,data=ads)
round(pairwise.t.test(clicks,type,data=ads),3)
library(multcomp)
install.packages("multcompView")
?multcompView
??multcompView
dif3 <- c(FALSE, FALSE, TRUE)
names(dif3) <- c("a-b", "a-c", "b-c")
multcompTs(dif3)
library(multcompView)
dif3 <- c(FALSE, FALSE, TRUE)
names(dif3) <- c("a-b", "a-c", "b-c")
multcompTs(dif3)
multcompLetters(dif3)
library(MASS)
multcompBoxplot(Postwt~Treat, data=anorexia)
setwd("~/Dropbox/courses/5210-2022/web-5210/psy5210/Projects/Chapter19")
ads.social <- 0:100/100
ads.normal <- 0:100/100
ads.animated <- 150:250/250
#Now, suppose that in an experiment, we have selected four advertisements of each type, and created a web page including all eight
#advertisements. In the study, each of 50 people either does or does not follow each link. Ignoring for the moment
#within-subject effects, we have a data set like this:
set.seed(100)
social.eff <- sample(ads.social,4)
normal.eff <- sample(ads.normal,4)
animated.eff <- sample(ads.animated,4)
##100 of each ad type
clicks <- c(runif(50*4)<social.eff, runif(50*4)<normal.eff, runif(50*4) <animated.eff) ##600 long; in
type <- as.factor(rep(c("A","B","C"),each=200))                                       ## 600 long
ad <- as.factor(rep(rep(1:4,50),3))
ad2 <- as.factor(c(rep(1:4,50),rep(11:14,50), rep(21:24,50)))                           ##ad id; in blocks of type.
##
sub <- as.factor(rep(rep(1:50,each=4),3))
contrasts(type) <- contr.sum(3)
contrasts(ad) <- contr.sum(length(unique(ad)))
contrasts(ad2) <- contr.sum(length(unique(ad2)))
contrasts(sub) <- contr.sum(50)
ads <- data.frame(sub=sub,type=type,ad=ad,ad2=ad2,clicks=clicks)
ads[1:10,]
##model without considering any random factors:
aggregate(ads$clicks,list(ads$type),mean)
##treat subject as a randomized factor:
model1 <-   aov(clicks~type+Error(sub/type),data=ads)
summary(model1)
##but we can't do a post-hoc test
pairwise.t.test(clicks,type,data=ads)
#type is a between-ad variable
model2 <- aov(clicks ~ type + Error(ad))
summary(model2)
#But what do you do?  According to Clark, you should use minF:
#fmin = F1 x F2/(F1+F2)
76.6*7.7/(76.6+7.7)
###########################################################
##Modern approach: mixed effects models:
library(nlme)
am2 <- lme(clicks ~ type, random = (~1|sub), data=ads)  #subject and randomized factor
am1 <- lme(clicks ~ 1, random = (~1|sub), data=ads)  #subject and randomized factor
library(lme4)
lmer2 <- lmer(clicks~type + (1|sub) ,data=ads)
summary(am2)
summary(lmer2)
library(multcomp)
contrasts(ads$type)
##The pairs you pick depend on how you did the contrasts!
## we used sum-to-zero contrasts. A is represented as 1 0, B is represented as 0 1,
## and C is represented as -1 -1.  So A versus B is row 1 minus Row 2;
## A versus C is row 1 plus row 3 plus row B, and  vs C is
#> contrasts(type)
#[,1] [,2]
#A    1    0
#B    0    1
#C   -1   -1
##contrasts are set up so that
contrasts <- rbind(
"Grand mean/intercept"=c(1,0,0),
"A alone"= c(1,1,0),
"B alone"= c(1,0,1),
"C alone"= c(1,-1,-1),
"A to B" = c(0,1,-1),
"A to C" = c(0,2,1),
"B to C" = c(0,1,2))
summary(glht(am2,contrasts))
summary(glht(am2,linfct=mcp(type="Tukey")))
contrasts
contrasts[5:7,]
summary(glht(am2,contrasts[5:7,]))
summary(glht(am2,linfct=mcp(type="Tukey")))
tukey <- glht(am2,
linfct=mcp(type="Tukey"))
print(tukey)
summary(tukey)
##compare to:
summary(glht(lmer2,contrast[5:7,]))
##compare to:
summary(glht(am2,contrast[5:7,]))
##compare to:
summary(glht(am2,contrasts[5:7,]))
##compare to:
summary(glht(am2,contrasts[5:7,]))
##compare to:
summary(glht(am2,contrasts[5:7,]))
##compare to:
summary(glht(am2,contrasts[5:7,]))
##compare to:
summary(glht(am2,contrasts[5:7,]))
##compare to:
summary(glht(am2,contrasts[5:7,]))
contrasts <- rbind(
"Grand mean/intercept"=c(1,0,0),
"A alone"= c(1,1,0),
"B alone"= c(1,0,1),
"C alone"= c(1,-1,-1),
"A to B" = c(0,1,-1),
"A to C" = c(0,2,1),
"B to C" = c(0,1,2))
contrasts
glht(am2,contrasts)
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
summary(glht(am2,contrasts))
am4 <- lme(clicks~type,random=list((~1|ad),(~1|sub)),data=ads) ##both
lmer4 <- lmer(clicks~type + (1|ad) + (1|sub),data=ads)
contrasts <- rbind(
"A to B" = c(0,1,-1),
"A to C" = c(0,2,1),
"B to C" = c(0,1,2))
summary(glht(lmer2,contrasts))
am4 <- lme(clicks~type,random=list((~1|ad),(~1|sub)),data=ads) ##both
lmer4 <- lmer(clicks~type + (1|ad) + (1|sub),data=ads)
summary(am4)
summary(lmer4)
summary(glht(lmer4,contrasts))
summary(glht(am5,contrasts))
?mcp
mcp
summary(glht(lmer4,contrasts))
summary(glht(lmer4, linfct=mcp(type="Tukey")))
cw <- ChickWeight
cw$logwt <- log(ChickWeight$weight)
contrasts(cw$Diet) <- contr.poly(levels(cw$Diet))
cw
lmer.cw0 <- lmer(logwt~Time*Diet + (1|Chick), data=cw )
summary(lmer.cw0)
Anova(lmer.cw0)
library(car)
Anova(lmer.cw0)
lmer.cw1 <- lmer(logwt~Time+Diet + (1|Chick),data=cw)
summary(lmer.cw1)
anova(lmer.cw0,lmer.cw1)
coef(lmer.cw0)
lmer.cw0
summary(lmer.cw0)
cor(coef(lmer.cw0))
coef(lmer.cw0)
coef(lmer.cw0)$Chick
cor(coef(lmer.cw0)$Chick)
cor(coef(lmer.cw0)$Chick,use="pairwise.complete")
lmer.cw2 <- lmer(logwt~Time*Diet + (1+Time|Chick),data=cw)
lmer.cw2b <- lmer(logwt~Time*Diet + (Time|Chick)+  (1|Chick) ,data=cw)
summary(glht(lmer.cw2, linfct=mcp(Diet="Tukey")))
summary(glht(lmer.cw2b, linfct=mcp(Diet="Tukey")))
summary(lmer.cw2)
Anova(lmer.cw2)
summary(lmer.cw2)
pr.cw <- profile(lmer.cw2)
confint(pr.cw)
pr.wc
pr.cw
pr.cw <- profile(lmer.cw2)
confint(pr.cw)
ads$clicks1 <- ads$clicks+0
library(ez)
ezm <- ezMixed(data=ads,dv=.(clicks1),
random=.(sub,ad),
fixed=.(type))
print(ezm)
mz <- (ezm$models$type$unrestricted)
summary(mz)
summary(glht(mz, linfct=mcp(type="Tukey")))
ezMixed
?ezMixed