R| 混合效应模型，lme4

## ref: A very basic tutorial for performing linear mixed effects analyses
## 分类变量
pitch = c(233,204,242,130,112,142)
sex = c(rep("female",3),rep("male",3))
my.df = data.frame(sex,pitch)
my.df
xmdl = lm(pitch ~ sex, my.df)
summary(xmdl)
mean(my.df[my.df$sex=="female",]$pitch)
## 连续变量
age = c(14,23,35,48,52,67)
pitch = c(252,244,240,233,212,204)
my.df = data.frame(age,pitch)
xmdl = lm(pitch ~ age, my.df)
summary(xmdl)

my.df$age.c = my.df$age - mean(my.df$age)
xmdl = lm(pitch ~ age.c, my.df)
summary(xmdl)

plot(fitted(xmdl),residuals(xmdl))
hist(residuals(xmdl))
qqnorm(residuals(xmdl))
dfbeta(xmdl)

###########################  混合效应模型  #####################################

library(lme4)
politeness=read.csv("/Users/Gavin/Downloads/politeness_data.csv")
# head(politeness), tail(politeness),summary(politeness), str(politeness), colnames(politeness)
which(is.na(politeness)==T)
boxplot(frequency ~ attitude*gender,
        col=c("white","lightgray"),politeness)
politeness.model = lmer(frequency ~ attitude +
                          (1|subject) + (1|scenario), data=politeness)
politeness.model

politeness.model = lmer(frequency ~ attitude +
                          gender + (1|subject) +
                          (1|scenario), data=politeness)
summary(politeness.model)
politeness.null = lmer(frequency ~ gender +
                         (1|subject) + (1|scenario), data=politeness,
                       REML=FALSE)
politeness.model = lmer(frequency ~ attitude +
                          gender +
                          (1|subject) + (1|scenario), data=politeness,
                        REML=FALSE)
anova(politeness.null,politeness.model)


coef(politeness.model)
politeness.model = lmer(frequency ~ attitude +
                          gender + (1+attitude|subject) +
                          (1+attitude|scenario),
                        data=politeness,
                        REML=FALSE)
coef(politeness.model)

politeness.null = lmer(frequency ~ gender +
                         (1+attitude|subject) + (1+attitude|scenario),
                       data=politeness, REML=FALSE)
# the likelihood ratio test
anova(politeness.null,politeness.model)





###############################################################################################
##       ref:https://ourcodingclub.github.io/2017/03/15/mixed-models.html#one              ##