# This code was developed to test mixed models for MUIR Wood zip paper


rm(list = ls())

library(ggplot2)
library(AICcmodavg)
library(lme4)
library(nlme)
library(AICcmodavg) # model averaging, perfered over MuMin
library(rcompanion) # for Tukey transformation of data
library(MuMIn)
library(ggeffects)

options(na.action = "na.fail")
# modified by CDF in Feb. 2021

#-------------------------------------------------------------------------------------------------------
# LOAD data
myFile <- "C:/Users/laf1021/Desktop/R Muir Zip/MUWOZipPaper20200622.csv" #choose.files() 
infile <- basename(myFile)
inpath <- dirname(myFile)
# READ IN DATA
dat <- read.csv(myFile)
#-------------------------------------------------------------------------------------------------------
#clean up missing data- 999
dat[dat==999] = NA
#get rid of rows with out all of the data 
dat = dat[complete.cases(dat$MEAN),]
dat = dat[complete.cases(dat$VisitorCount),]
#------------------------------------------------------------------------------------------------------

#-------------------------------------------------------------------------------------------------------

#-------------------------------------------------------------------------------------------------------

dat$Pleasantne = as.numeric( dat$Pleasantne )  
hist(dat$Pleasantne)



### I tried all kinds of transformations for left or negative-skew data. Below seems to achieve the best for improving the residuals of the model
### this Tukey transformation does make the response variable normally distributed, but it does improve residuals of the models below
dat$PleasT <- transformTukey(dat$Pleasantne)

hist(dat$PleasT)
#-------------------------------------------------------------------------------------------------------



# Distribution of the response variables...
par(mfcol=c(2,2), mar=c(5,4,2,0.5), oma=c(0.5,0.5,0.5,0.5))
#all clips model
hist(dat$Pleasantne,main = "Pleasantess",xlab = "", freq = FALSE)

#-------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------
# characterize variables correctly... as.numeric, as.factor,
dat$MEAN = as.numeric( dat$MEAN)
dat$HourlyLEQ = as.numeric ( dat$HourlyLEQ)
dat$HourlyL50 = as.numeric ( dat$HourlyL50)
dat$AbleHear = as.numeric( dat$AbleHear)
dat$NSS_SF = as.numeric ( dat$NSS_SF)
dat$RC_SoundMotive = as.numeric ( dat$RC_SoundMotive)
dat$VisitorCount = as.numeric( dat$VisitorCount)
dat$DummyQuiet = as.numeric( dat$DummyQuiet)
dat$Time = as.numeric( dat$Time)
dat$date = as.numeric( dat$date)
dat$Day = as.numeric( dat$Day)

## make time a factor for consideration as a random intercept
dat$Time <- as.factor(dat$Time)


#-------------------------------------------------------------------------------------------------------
# BUILD each model 
library(spaMM)

## use marginal AIC to rank models

global <- fitme(Pleasantne ~ MEAN+NSS_SF+AbleHear+HourlyL50+VisitorCount + Matern(1|Xcoord + xcoor) + AR1(1|Day), data = dat,fixed=list(nu=0.5))
summary(global)
get_any_IC(global) # 1003.19


globalnoTime <- fitme(Pleasantne ~ MEAN+NSS_SF+AbleHear+HourlyL50+VisitorCount + Matern(1|Xcoord + xcoor), data = dat,fixed=list(nu=0.5))
summary(globalnoTime)
get_any_IC(globalnoTime) # 999.32

globalnoSpace <- fitme(Pleasantne ~ MEAN+NSS_SF+AbleHear+HourlyL50+VisitorCount + AR1(1|Day), data = dat)
summary(globalnoSpace)
get_any_IC(globalnoSpace) # 999.186

globalnoSpaceTime <- fitme(Pleasantne ~ MEAN+NSS_SF+AbleHear+HourlyL50+VisitorCount, data = dat)
summary(globalnoSpaceTime)
get_any_IC(globalnoSpaceTime) # 995.32

globalnoSpaceTimeR <- fitme(Pleasantne ~ MEAN+NSS_SF+AbleHear+HourlyL50+VisitorCount + (1|Day/Time), data = dat)
summary(globalnoSpaceTimeR)
get_any_IC(globalnoSpaceTimeR) ## 999.2







#-------------------------------------------------------------------------------------------------------
# switch to linear models because no evidence to need correlation structure for space or time OR random effects


# Global model
m7 <-	lm(Pleasantne ~ MEAN + NSS_SF + AbleHear+VisitorCount+HourlyL50, data = dat)
m7s <- lm(PleasT ~ MEAN + NSS_SF + AbleHear+VisitorCount+HourlyL50, data = dat)

# nontransformed version
drm7 <- dredge(m7)
drm7

# transformed version
drm7s <- dredge(m7s)
drm7s

#Build Linear models and use Stepwise/AiC 3/transformed pleasantness
Cand.models <- list()
#hypothesis: 

Cand.models[[1]] <-	lm(PleasT ~ MEAN,  data = dat)
summary(Cand.models[[1]])
AIC(Cand.models[[1]])
Cand.models[[2]] <-	lm(PleasT ~ MEAN+NSS_SF,   data = dat)
summary(Cand.models[[2]])
AIC(Cand.models[[2]])
Cand.models[[3]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear,   data = dat)
summary(Cand.models[[3]])
AIC(Cand.models[[3]])
Cand.models[[4]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive,  data = dat)
summary(Cand.models[[4]])
AIC(Cand.models[[4]])
Cand.models[[5]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive+DummyQuiet,  data = dat)
summary(Cand.models[[5]])
AIC(Cand.models[[5]])
Cand.models[[6]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive+VisitorCount,  data = dat)
summary(Cand.models[[6]])
AIC(Cand.models[[6]])
Cand.models[[7]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive+HourlyL50,  data = dat)
summary(Cand.models[[7]])
AIC(Cand.models[[7]])
Cand.models[[8]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive+AbleHear:DummyQuiet,  data = dat)
summary(Cand.models[[8]])
AIC(Cand.models[[8]])
Cand.models[[9]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive+AbleHear:HourlyL50,  data = dat)
summary(Cand.models[[9]])
AIC(Cand.models[[9]])
Cand.models[[10]] <-	lm(PleasT ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive+AbleHear:VisitorCount,  data = dat)
summary(Cand.models[[10]])
AIC(Cand.models[[10]])


#Best Model is 4!!
summary(Cand.models[[4]])
AIC(Cand.models[[4]])

# top models
#non-transformed
m3 <-	lm(Pleasantne ~ MEAN+NSS_SF+AbleHear+RC_SoundMotive, data = dat)
summary(m3)
## plot diagnostics
qqnorm(m3$residuals) ### this does not look great
qqline(m3$residuals)
hist(m3$residuals)


m3s <-	lm(PleasT ~ MEAN + NSS_SF + AbleHear + RC_SoundMotive, data = dat)
summary(m3s)

qqnorm(m3s$residuals) ### this looks better than non-transformed
qqline(m3s$residuals)
hist(m3s$residuals) # also looks better

par(mfrow=c(2,2))


# ## plots of model 3
DFsm3c <- ggpredict(m3, terms = c("MEAN"), back.transform =T, type ="fe")

C1 <- plot(DFsm3c, colors="eight", add.data=T)  + labs(x="Sound level at home zip (dB)", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black")) +
  labs(subtitle = "(c)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
C1


DFsm3a <- ggpredict(m3, terms = c("NSS_SF"), back.transform =T, type ="fe")

A1 <- plot(DFsm3a, colors="eight", add.data=T)  + labs(x="Noise sensitivity Index", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))+
  labs(subtitle = "(a)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
A1

DFsm3b <- ggpredict(m3, terms = c("AbleHear"), back.transform =T, type ="fe")

B1 <- plot(DFsm3b, colors="eight", add.data=T)  + labs(x="Noise interference Index", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))+
  labs(subtitle = "(b)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
B1



library("gridExtra")
grid.arrange(A1, B1, C1,  
             ncol = 3, nrow = 1)
g <- arrangeGrob(A1, B1, C1,  
                 ncol = 3, nrow = 1)

ggsave("MarginalEffects.pdf",g,width = 20, height = 7.5, units = "cm")



#### plots of m3s



# ## plots of model 3, note that these are the transformed versions of pleasantness. 
DFsm3sc <- ggpredict(m3s, terms = c("MEAN"), back.transform =T, type ="fe")

C1s <- plot(DFsm3sc, colors="eight", add.data=T)  + labs(x="Neighborhood sound level (dB)", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black")) +
  labs(subtitle = "(c)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
C1s


DFsm3sa <- ggpredict(m3s, terms = c("NSS_SF"), back.transform =T, type ="fe")

A1s <- plot(DFsm3sa, colors="eight", add.data=T)  + labs(x="Noise sensitivity Index", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))+
  labs(subtitle = "(a)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
A1s

DFsm3sb <- ggpredict(m3s, terms = c("AbleHear"), back.transform =T, type ="fe")

B1s <- plot(DFsm3sb, colors="eight", add.data=T)  + labs(x="Noise interference Index", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))+
  labs(subtitle = "(b)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
B1s

DFsm3sb <- ggpredict(m3s, terms = c("RC_SoundMotive"), back.transform =T, type ="fe")

D1s <- plot(DFsm3sb, colors="eight", add.data=T)  + labs(x="Sound Motivation Index", y ="Pleasantness", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))+
  labs(subtitle = "(d)")+ theme(plot.subtitle = element_text(face="bold", size = 14))
D1s

library("gridExtra")
grid.arrange(A1s, B1s, C1s, D1s,  
             ncol = 4, nrow = 1)
g <- arrangeGrob(A1s, B1s, C1s,  
                 ncol = 4, nrow = 1)

ggsave("MarginalEffectsTransformed.pdf",g,width = 20, height = 7.5, units = "cm")

############################The performance package checks for multicollinearity and the function is check_collinearity(model).##
#### how is MEAN, NSS_SF and AbleHear related?
## visitor count has a small negative influence on noise interference
mAH <-	lm(AbleHear ~ MEAN + NSS_SF + RC_SoundMotive, data = dat)
summary(mAH)

#Checked for Interactions between AbleHear and other variables you suggested on 3.7.23 - looks like very small/positive 
mAH2 <-	lm(AbleHear ~ HourlyL50 + DummyQuiet + VisitorCount, data = dat)
summary(mAH2)

mNS <-	lm(NSS_SF ~ MEAN + AbleHear + RC_SoundMotive, data = dat)
summary(mNS)

qqnorm(resid(mNS))
qqline(resid(mNS))
hist(resid(mNS)) ## these all look good


DFmNS <- ggpredict(mNS, terms = c("MEAN"), back.transform =T, type ="fe")

NS1 <- plot(DFmNS, colors="eight", add.data=T)  + labs(x="Neighborhood sound level (dB)", y ="Noise sensitivity", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))
NS1
#### how is MEAN, NSS_SF and AbleHear related?
## visitor count has a small negative influence on noise interference
mAH <-	lm(AbleHear ~ MEAN + NSS_SF + RC_SoundMotive, data = dat)
summary(mAH)


mNS <-	lm(NSS_SF ~ MEAN + AbleHear + RC_SoundMotive, data = dat)
summary(mNS)

qqnorm(resid(mNS))
qqline(resid(mNS))
hist(resid(mNS)) ## these all look good


DFmNS <- ggpredict(mNS, terms = c("MEAN"), back.transform =T, type ="fe")

NS1 <- plot(DFmNS, colors="eight", add.data=T)  + labs(x="Nieghborhood sound level (dB)", y ="Noise sensitivity", title ="") + theme_classic() + theme(legend.position = "none") + theme(panel.grid.minor = element_blank())+ theme(axis.line = element_line(colour = "black"))
NS1

#Check for collinearity 
library(performance)

check_collinearity(Cand.models[[4]])