##################
#  HLM Models    #
##################

#This code generates HLM models and spotlights for Studies 1a-1c and 2a-2c

library(tidyverse)
library(Hmisc)
library(lme4)
library(lmerTest)
library(sjPlot)


#### Inputs ####
load("AllData.Rdata")
StudyList <- c("AgeProb","AgeMiss","AgeQuant","FoodProb","FoodMiss","FoodQuant")


#### Remove outliers ####
data <- AllData %>%
  filter(OC_Change <5000 & OC_Change > -5000 & ExpImprove <5000 & ExpImprove >-5000) %>% 
  group_by(Study) %>%
  mutate_at(vars(OC_Change, ExpImprove, Improve, DistIA), list(scaled = scale)) %>% 
  ungroup() 
AllData5sd <- data %>% 
  filter(abs(OC_Change_scaled)<5 & abs(ExpImprove_scaled)<5 & abs(Improve_scaled)<5 & abs(DistIA_scaled)<5)


##### Table 4. HLM regression table #####
#All IVs are mean centered (indicated by "c" at end of variable name)
hlmOCcentered <- function (StimElicit) {
  d <- filter(AllData5sd,Study == StimElicit)
  Lev1Fixed <- lmer(OC_Change ~ Truthc + DistIAc + Rankc + DistIAc*Rankc 
                    + I(DistIAc^2)  
                    + I(Rankc^2)  
                    + (1 |Subject)+ (1|StimNum), d) 
}
m1<-hlmOCcentered("AgeProb")
m2<-hlmOCcentered("AgeMiss")
m3<-hlmOCcentered("AgeQuant")
m4<-hlmOCcentered("FoodProb")
m5<-hlmOCcentered("FoodMiss")
m6<-hlmOCcentered("FoodQuant")
tab_model(m1,m2,m3,m4,m5,m6, show.se = TRUE, show.ci = FALSE, digits = 5, show.aic = TRUE, show.loglik = TRUE, show.obs = TRUE, dv.labels = StudyList)


#Model comparisons
#This compares two different models
mcfn <- function (regmodel1, regmodel2) {
      for (i in 1:6) {
          StimElicit <- StudyList[i]
          cat("\n",
                        "--------------------------------------------------", "\n", 
                        StimElicit, "\n",
                        "--------------------------------------------------", "\n")
          print(anova(regmodel1(StimElicit),regmodel2(StimElicit)))
        }
    }    
hlmOCNoRankNoDist <- function (StimElicit) {
    d <- filter(AllData5sd,Study == StimElicit)
    Lev1Fixed <- lmer(OC_Change ~ Truthc +  (1 |Subject)+ (1|StimNum), d) 
  }  
hlmOCBothQuad <- function (StimElicit) {
    d <- filter(AllData5sd,Study == StimElicit)
    Lev1FixedQuad <- lmer(OC_Change ~ Truthc + DistIAc + Rankc + DistIAc*Rankc 
                                                   + I(DistIAc^2) + I(Rankc^2) 
                                                   + (1 |Subject)+ (1|StimNum), d) 
    }
mcfn(hlmOCBothQuad,hlmOCNoRankNoDist)



##### Table 5. Spotlights #####
#Note: This is series of nested functions that:
#1. Select the data for the relevant study and spotlight values
#2. Assign a functional form for the regression
#3. Aggregate the intercept coefficient estimates into a table


#1. Select data
selectdata <- function (DataSet, StimElicit, DV, Dist, AdvQuality){
  d <- filter(DataSet, Study == StimElicit)%>%
    rowwise() %>%
    mutate(
      DepVar = ifelse(DV == "OC_Change", OC_Change,
                      ifelse(DV == "ExpImprove", ExpImprove,
                             ifelse(DV == "Improve", Improve, NA))),
      DistVar = ifelse(Dist == 0, DistIA,
                       ifelse(Dist == 10, Dist10c,
                              ifelse(Dist == 50, Dist50c,
                                     ifelse (Dist == 90, Dist90c, NA)))),
      RankVar = ifelse(AdvQuality == "High", Rank25c,
                       ifelse(AdvQuality == "Mean",Rankc,
                              ifelse(AdvQuality == "Low", Rank175c, NA)))
    ) %>%
    ungroup()
}


#2. Assign functional form
#BothQuad
bothquad <- function (d){
  Lev1Fixed <- lmer(DepVar ~ Truthc + DistVar + RankVar + DistVar*RankVar 
                    + I(DistVar^2) + I(RankVar^2) 
                    + (1 |Subject)+ (1|StimNum), d) 
}

#3. Aggregate intercept coefficients
#Inputs to the function
StudyList <- c("AgeProb","AgeMiss","AgeQuant","FoodProb","FoodMiss","FoodQuant")
DVList <- c("OC_Change", "ExpImprove", "Improve")
DistList <- c(10, 90)
AdvQualityList <- c("High", "Low")
DVOrder <- 1:3
DVDetails <-data.frame (DVList, DVOrder)
EstList <- c("Estimate", "Std. Error","df", "t value","Pr(>|t|)")
EstOrder <- 1:5
EstDetails <- data.frame (EstList, EstOrder)

cofn <- function (DataSet, regmodel, modelname) {
  codf <- data.frame(matrix(ncol=0,nrow=0)) #empty data frame for appending
  for (i in 1:6) {
    StimElicit <- StudyList[i]
    for (j in 1:3) {
      DV <- DVList[j]
      for (k in 1:2) {
        Dist <- DistList[k]
        for (l in 1:2) {
          AdvQuality <- AdvQualityList[l]
          co <- as.data.frame(summary(regmodel(selectdata(DataSet, StimElicit, DV, Dist, AdvQuality)))$coefficients)[1,] %>% 
            mutate(Study = StimElicit,
                   DV = DV,
                   Dist = Dist,
                   AdvQuality = AdvQuality)
          codf<- bind_rows(codf,co)
        }
      }
    }
  } 
  codf <-as.data.frame(codf) %>% 
    gather(key = "var", value = "val", -c(Study, DV, Dist, AdvQuality)) %>%
    spread(key = "Study", value = "val") %>%  
    inner_join(DVDetails, by = c("DV" = "DVList")) %>%
    inner_join(EstDetails, by = c("var" = "EstList")) %>%
    arrange(AdvQuality, Dist, DVOrder, EstOrder) %>% 
    select(AdvQuality, Dist, DV, var, AgeProb,AgeMiss,AgeQuant,FoodProb,FoodMiss,FoodQuant)
}



#Produce df of all intercept coefficients
AllCoeffsBothQuadNoOut5sd <- cofn(AllData5sd,bothquad,"bothquad")
write_csv(AllCoeffsBothQuadNoOut5sd, "AllCoeffsBothQuadNoOut5sd.csv")

