library(ggplot2)       # for 数据可视化
library(lme4)          # for 混合效应模型（多层模型）的构建
library(lmerTest)      # for 计算P值
library(psych)         # for 描述统计
library(plyr)          # for 数据整理
library(effects)       # for 检验调节效应
library(interactions)  # for 调节效应可视化

pdata = read.csv("person-data.csv")
ddata = read.csv("daily-data.csv")
pdata <- pdata[ ,c("id","bfi_n")]
ddata <- ddata[ ,c("id","day","negaff","pss")]
head(pdata)

head(ddata)

ddata$stress <- 4 - ddata$pss

psych::describe(ddata$stress)

ggplot(data=ddata, aes(x=stress)) +
  geom_histogram(fill="white", color="black",bins=19) +
  labs(x = "压力")

# 计算个体均值
# negaff_trait虽然不是模型中的变量， 但是我们后期可能会用到它，所以也在这里一并计算了
personmeans <- ddply(ddata, "id", summarize,
                       stress_trait = mean(stress, na.rm=TRUE),
                       negaff_trait = mean(negaff, na.rm=TRUE))
head(personmeans)
describe(personmeans)

# 将新生成的个体特质变量合并到pdata数据框, pdata 是被试样本
pdata <- merge(pdata, personmeans, by="id")                                              
# 后面我们会分析调节效应， 对变量中心化是必然的步骤，这一步就是变量中心化
pdata$bfi_n_c <- scale(pdata$bfi_n,center=TRUE,scale=FALSE)
pdata$stress_trait_c <- scale(pdata$stress_trait,center=TRUE,scale=FALSE)
# 个体数据的描述性统计
describe(pdata)

ddata_long <- merge(ddata,pdata,by="id")

# 计算状态变量，即减去个体均值
ddata_long$stress_state <- ddata_long$stress - ddata_long$stress_trait
ddata_long$negaff_state <- ddata_long$negaff - ddata_long$negaff_trait

# 查看数据
head(ddata_long)

#faceted plot
ggplot(data=ddata_long[which(ddata_long$id <= 125),], aes(x=stress_state,y=negaff)) +
  geom_point() +
  stat_smooth(method="lm", fullrange=TRUE) +
  xlab("Stress State") + ylab("Negative Affect (Continuous)") + 
  facet_wrap( ~ id) +
  theme(axis.title=element_text(size=16),
        axis.text=element_text(size=14),
        strip.text=element_text(size=14))

#unconditional means model
model0_fit <- lmer(formula = negaff ~ 1 + (1|id), 
              data=ddata_long,
              na.action=na.exclude)
summary(model0_fit)

# 我们使用 VarCorr 函数提取方差 ， id这一行对应的是被试内方差， Residual对应的是残差方差（不可以被被试差异所解释的方差
VarCorr(model0_fit)

randEffs <- as.data.frame(VarCorr(model0_fit))
randEffs

u0v = randEffs[1,4]
ev = randEffs[2,4]

icc <- u0v / (u0v+ev)
icc

model1 <- lmer(formula = negaff ~ 1 + day + stress_trait_c + 
                      stress_state + stress_state:stress_trait_c + 
                      (1 + stress_state|id), 
                    data=ddata_long,
                    na.action=na.exclude)
summary(model1)

# 保存模型的预测结果
ddata_long$pred_m1 <- predict(model1)
head(ddata_long)

confint(model1)

ggplot(data=personmeans, aes(x=stress_trait, y=negaff_trait, group=factor(id)), legend=FALSE) +
  geom_point(colour="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
  xlab("Trait Stress") + ylab("Trait Negative Affect") +
  theme_classic() +
  theme(axis.title=element_text(size=16),
        axis.text=element_text(size=12),
        plot.title=element_text(size=16, hjust=.5)) +
  ggtitle("Between-Person Association Plot\nTrait Stress & Negative Affect")

ggplot(data=ddata_long, aes(x=stress_state, y=negaff_state, group=factor(id), colour="gray"), legend=FALSE) +
  geom_smooth(method=lm, se=FALSE, fullrange=FALSE, lty=1, size=.5, color="gray40") +
  geom_smooth(aes(group=1), method=lm, se=FALSE, fullrange=FALSE, lty=1, size=2, color="blue") +
  xlab("Stress State") + ylab("Predicted State Negative Affect") +
  theme_classic() +
  theme(axis.title=element_text(size=18),
        axis.text=element_text(size=14),
        plot.title=element_text(size=18, hjust=.5)) +
  ggtitle("Within-Person Association Plot\nPerceived Stress & Negative Affect")

# fit model
model2 <- lmer(formula = negaff ~ 1 + day + stress_trait_c + 
                      bfi_n_c + stress_trait_c:bfi_n_c +
                      stress_state + stress_state:stress_trait_c + 
                      stress_state:bfi_n_c + stress_state:stress_trait_c:bfi_n_c + 
                      (1 + stress_state|id),
                    data=ddata_long,
                    na.action=na.exclude)
#Look at results
summary(model2)

ddata_long$pred_m2 <- predict(model2)
head(ddata_long)

describe(ddata_long$bfi_n_c)
describe(ddata_long$stress_state)

#calculate effect
effects_model2 <- effect(term="bfi_n_c*stress_state", mod=model2,xlevels=list(bfi_n_c=c(-0.96, +0.96), stress_state=c(-0.49,+0.49)))
summary(effects_model2)

#convert to dataframe
effectsdata <- as.data.frame(effects_model2)
#plotting the effect evaluation (with standard error ribbon)
ggplot(data=effectsdata, aes(x=stress_state, y=fit, group=bfi_n_c), legend=FALSE) + 
  geom_point() +
  geom_line() +
  #geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.3) +
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.15) +
  xlab("Stress State") + xlim(-2,2) +
  ylab("Predicted Negative Affect") + ylim(1,7) +
  ggtitle("Differences in Stress Reactivity across Neuroticism")

johnson_neyman(model=model2, pred=stress_state, modx=bfi_n_c)

BIC(logLik(model2))

logLik(logLik(model2))

BIC(logLik(model2))