dat<- read.csv(file='AFC Advertisement Final Data.csv',header=TRUE,dec=',',row.names = NULL,sep=';')


datc <- dat[complete.cases(dat),]

library(lme4)

library(gamm4)

library(itsadug)

library(mgcv)

library(ggplot2)

library(caret)
options(max.print = 100000)


nearZeroVar(dat, names = TRUE)


#Models for Table 2: 
#Prediction of Ad Likeability (Mod1A), 
#Brand Likeability Change (Mod1B), 
#and Purchase Intention Change (Mod1C)
#based on relevant self-report ratings

Mod1A <- gamm4(ZAdLike_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus) 
               + s(Zhappy_re)
               
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod1B <- gamm4(ZBrandLikeDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(ZAdLike_re)
               + s(Zhappy_re)
               
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod1C <- gamm4(ZBrandPIDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)  
               + s(ZBrandLikeDiff)
               + s(ZAdLike_re)
               + s(Zhappy_re)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

#Models for Table 3: 
#Prediction of joy ratings based on mean (Mod2A) and peak AU activity (Mod3A)


Mod2A <- gamm4(Zhappy_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod3A <- gamm4(Zhappy_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

#Models for Table 4: 
#Prediction of advertisement likeability ratings based on 
#mean AU Activity without self-report ratings (Mod2B)
#mean AU Activity with self-report ratings (Mod4B)
#peak AU Activity without self-report ratings (Mod3B)
#peak AU Activity with self-report ratings (Mod5B)


#mean and peak AU Activity with and without self-report ratings


Mod2B <- gamm4(ZAdLike_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               
               ,random = ~ (1|Participant)  + (1|Experiment)
               ,data=dat)

Mod4B <- gamm4(ZAdLike_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus) 
               + s(Zhappy_re)
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod3B <- gamm4(ZAdLike_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod5B <- gamm4(ZAdLike_re ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus) 
               + s(Zhappy_re)
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)


#Models for Table 5: 
#Prediction of brand likeability change ratings based on 
#mean AU Activity without self-report ratings (Mod2C)
#mean AU Activity with self-report ratings (Mod8C)
#peak AU Activity without self-report ratings (Mod3C)
#peak AU Activity with self-report ratings (Mod7C)

Mod2C <- gamm4(ZBrandLikeDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod8C <- gamm4(ZBrandLikeDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(ZAdLike_re)
               + s(Zhappy_re)
               
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod3C <- gamm4(ZBrandLikeDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod7C <- gamm4(ZBrandLikeDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(ZAdLike_re)
               + s(Zhappy_re)
               
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)


#Models for Table 6: 
#Prediction of brand likeability change ratings based on 
#mean AU Activity without self-report ratings (Mod2D)
#mean AU Activity with self-report ratings (Mod8D)
#peak AU Activity without self-report ratings (Mod3D)
#peak AU Activity with self-report ratings (Mod7D)

Mod2D <- gamm4(ZBrandPIDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod8D <- gamm4(ZBrandPIDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(ZAdLike_re)
               + s(Zhappy_re)
               + s(ZBrandLikeDiff)
               
               + s(Mean_AU04BrowLowerer_A)
               + s(Mean_AU05UpperLidRaiser_A)
               + s(Mean_AU06CheekRaiser_A)
               + s(Mean_AU07LidTightener_A)
               + s(Mean_AU12_A)
               + s(Mean_AU14Dimpler_A)
               + s(Mean_AU15LipCornerDepressor_A)
               + s(Mean_AU17ChinRaiser_A)
               + s(Mean_AU23LipTightener_A)
               + s(Mean_AU24LipPressor_A)
               + s(Mean_AU25LipsPart_A)
               + s(Mean_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)


Mod3D <- gamm4(ZBrandPIDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)

Mod7D <- gamm4(ZBrandPIDiff ~ 
                 + as.factor(Gender)
               + as.factor(Stimulus)                 
               + s(ZAdLike_re)
               + s(Zhappy_re)
               + s(ZBrandLikeDiff)
               
               + s(Max_AU04BrowLowerer_A)
               + s(Max_AU05UpperLidRaiser_A)
               + s(Max_AU06CheekRaiser_A)
               + s(Max_AU07LidTightener_A)
               + s(Max_AU12_A)
               + s(Max_AU14Dimpler_A)
               + s(Max_AU15LipCornerDepressor_A)
               + s(Max_AU17ChinRaiser_A)
               + s(Max_AU23LipTightener_A)
               + s(Max_AU24LipPressor_A)
               + s(Max_AU25LipsPart_A)
               + s(Max_AU43EyesClosed_A)
               ,random = ~ (1|Participant) + (1|Experiment)
               ,data=dat)


#Summaries

#Models for Table 2: 
summary(Mod1A$gam) 
summary(Mod1B$gam)  
summary(Mod1C$gam) 

AIC(Mod1A$mer, Mod1B$mer, Mod1C$mer)
BIC(Mod1A$mer, Mod1B$mer, Mod1C$mer)

#Models for Table 3: 
summary(Mod2A$gam) 
summary(Mod3A$gam) 

AIC(Mod2A$mer, Mod3A$mer)
BIC(Mod2A$mer, Mod3A$mer)

#Models for Table 4: 
summary(Mod2B$gam)  
summary(Mod4B$gam)  
summary(Mod3B$gam)  
summary(Mod5B$gam)  

AIC(Mod2B$mer, Mod4B$mer, Mod3B$mer, Mod5B$mer)
BIC(Mod2B$mer, Mod4B$mer, Mod3B$mer, Mod5B$mer)

#Models for Table 5: 
summary(Mod2C$gam) 
summary(Mod8C$gam) 
summary(Mod3C$gam) 
summary(Mod7C$gam) 

AIC(Mod2C$mer, Mod8C$mer, Mod3C$mer, Mod7C$mer)
BIC(Mod2C$mer, Mod8C$mer, Mod3C$mer, Mod7C$mer)

#Models for Table 6: 
summary(Mod2D$gam) 
summary(Mod8D$gam) 
summary(Mod3D$gam) 
summary(Mod7D$gam) 

AIC(Mod2D$mer, Mod8D$mer, Mod3D$mer, Mod7D$mer)
BIC(Mod2D$mer, Mod8D$mer, Mod3D$mer, Mod7D$mer)

#############################
PLOTS
#width=1674, height=1572, 

tiff(filename="Figure 2.tiff", width=2500, height=2300, res=300)
par(mfrow=c(3,4))
plot(Mod2A$gam,select=5, ylim=c(-1,1), xlim=c(0,0.8), main="1A", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Mean", ylab="Joy Rating") 
abline(h=0)
plot(Mod2B$gam,select=5, ylim=c(-1,1), xlim=c(0,0.8), main="1B", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Mean", ylab="Ad Like")
abline(h=0)
plot(Mod2C$gam,select=5, ylim=c(-1,1), xlim=c(0,0.8), main="1C", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Mean", ylab="Brand Like Change")
abline(h=0)
plot(Mod2D$gam,select=5, ylim=c(-1,1), xlim=c(0,0.8), main="1D", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Mean", ylab="Purchase Intention Change")
abline(h=0)


plot(Mod2A$gam,select=3, ylim=c(-2,2), xlim=c(0,0.8), main="2A", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Mean", ylab="Joy Rating") 
abline(h=0)
plot(Mod2B$gam,select=3, ylim=c(-2,2), xlim=c(0,0.8), main="2B", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Mean", ylab="Ad Like")
abline(h=0)
plot(Mod2C$gam,select=3, ylim=c(-2,2), xlim=c(0,0.8), main="2C", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Mean", ylab="Brand Like Change")
abline(h=0)
plot(Mod2D$gam,select=3, ylim=c(-2,2), xlim=c(0,0.8), main="2D", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Mean", ylab="Purchase Intention Change")
abline(h=0)

plot(Mod2A$gam,select=6, ylim=c(-3,3), xlim=c(0,0.8), main="3A", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 14 Mean", ylab="Joy Rating") 
abline(h=0)
plot(Mod2B$gam,select=6, ylim=c(-3,3), xlim=c(0,0.8), main="3B", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 14 Mean", ylab="Ad Like")
abline(h=0)
plot(Mod2C$gam,select=6, ylim=c(-3,3), xlim=c(0,0.8), main="3C", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 14 Mean", ylab="Brand Like Change")
abline(h=0)
plot(Mod2D$gam,select=6, ylim=c(-3,3), xlim=c(0,0.8), main="3D", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 14 Mean", ylab="Purchase Intention Change")
abline(h=0)

dev.off()


tiff(filename="Figure 3.tiff", width=2500, height=2300, res=300)
par(mfrow=c(3,4))
plot(Mod3A$gam,select=5, ylim=c(-1,1), main="1A", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Peak", ylab="Joy Rating") 
abline(h=0)
plot(Mod3B$gam,select=5, ylim=c(-1,1), main="1B", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Peak", ylab="Ad Like")
abline(h=0)
plot(Mod3C$gam,select=5, ylim=c(-1,1), main="1C", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Peak", ylab="Brand Like Change")
abline(h=0)
plot(Mod3D$gam,select=5, ylim=c(-1,1), main="1D", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 12 Peak", ylab="Purchase Intention Change")
abline(h=0)


plot(Mod3A$gam,select=3, ylim=c(-1,1), main="2A", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Peak", ylab="Joy Rating") 
abline(h=0)
plot(Mod3B$gam,select=3, ylim=c(-1,1), main="2B", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Peak", ylab="Ad Like")
abline(h=0)
plot(Mod3C$gam,select=3, ylim=c(-1,1), main="2C", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Peak", ylab="Brand Like Change")
abline(h=0)
plot(Mod3D$gam,select=3, ylim=c(-1,1), main="2D", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 06 Peak", ylab="Purchase Intention Change")
abline(h=0)

plot(Mod3A$gam,select=1, ylim=c(-1,1), main="3A", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 04 Peak", ylab="Joy Rating") 
abline(h=0)
plot(Mod3B$gam,select=1, ylim=c(-1,1), main="3B", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 04 Peak", ylab="Ad Like")
abline(h=0)
plot(Mod3C$gam,select=1, ylim=c(-1,1), main="3C", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 04 Peak", ylab="Brand Like Change")
abline(h=0)
plot(Mod3D$gam,select=1, ylim=c(-1,1), main="3D", shade=TRUE, rug=TRUE, shade.col='grey80', xlab="Action Unit 04 Peak", ylab="Purchase Intention Change")
abline(h=0)

dev.off()