#Variance analysis
setwd ("C:/Users/lenovo/Desktop")
df<- read.csv("C:/Users/lenovo/Desktop/Data.csv")

##### T-test#######
t.test(df$DE~df$Residence)
t.test(df$IDE~df$Residence)
t.test(df$VE~df$Residence)
t.test(df$CC~df$Residence)
t.test(df$EC~df$Residence)
t.test(df$PE~df$Residence)
t.test(df$PN~df$Residence)

##### ANOVA #######
aov1 <- aov(DE~Location, df) 
aov1
summary(aov1)

aov2 <- aov(IDE~Location, df) 
aov2
summary(aov2)

aov3 <- aov(VE~Location, df) 
aov3
summary(aov3)

aov4 <- aov(CC~Location, df) 
aov4
summary(aov4)

aov5 <- aov(EC~Location, df) 
aov5
summary(aov5)

aov6 <- aov(PE~Location, df) 
aov6
summary(aov6)

aov7 <- aov(PN~Location, df) 
aov7
summary(aov7)
##### Tukey honestly significant difference (HSD) multiple comparisons #######
fit1 <- aov1
summary(fit1)
TukeyHSD(fit1)
plot(TukeyHSD(fit1))

fit2 <- aov2
summary(fit2)
TukeyHSD(fit2)
plot(TukeyHSD(fit2))

fit3 <- aov3
summary(fit3)
TukeyHSD(fit3)
plot(TukeyHSD(fit3))

fit4 <- aov4
summary(fit4)
TukeyHSD(fit4)
plot(TukeyHSD(fit4))

fit5 <- aov5
summary(fit5)
TukeyHSD(fit5)
plot(TukeyHSD(fit5))

fit6 <- aov6
summary(fit6)
TukeyHSD(fit6)
plot(TukeyHSD(fit6))

fit7 <- aov7
summary(fit7)
TukeyHSD(fit7)
plot(TukeyHSD(fit7))

####### The interaction effect########
###### Two-way ANOVA######
library(gplots)
library(HH) 

## 1.DE
interaction2wt(df$DE~df$Residence*df$Location)

## 2.IDE
interaction2wt(df$IDE~df$Residence*df$Location)

## 3.VE
interaction2wt(df$VE~df$Residence*df$Location)

## 4.CC
interaction2wt(df$CC~df$Residence*df$Location) 

## 5. EC
interaction2wt(df$EC~df$Residence*df$Location) 

## 6.PE
interaction2wt(df$PE~df$Residence*df$Location) 

##7.PN
interaction2wt(df$PN~df$Residence*df$Location) 


Mixed linear model

library(car)

##### 1.Pro-environmental behavior (PE)#######

fit<-lm(PE~ DE+IDE+VE+DE*IDE+DE*VE+IDE*VE++CC+EC+ Gender+Age+Residence+Location+DE*Residence+IDE*Residence+VE*Residence,data = df)
summary(fit)
vif(fit)
plot(fit)
fit_step <- step(fit) 
summary(fit_step)
vif(fit_step)

## Multiple Regression                                               
## Use the full model as a star

fit<-lm(PE~ DE+IDE+VE+DE*IDE+DE*VE+IDE*VE++CC+EC+ Gender+Age+Residence+Location+DE*Residence+IDE*Residence+VE*Residence,data = df);
summary (fit)

## Update model
fit1 <- step (fit) 
summary(fit1)

# Evaluate model
# multicollinearity

vif(lm(PE~DE+IDE+VE+CC+ EC+Gender+Age+Residence+Location , data = df))
final <- lm(PE ~ DE + IDE + VE + CC+ EC + Gender + Age + Residence + Location, data = df)

# Independence, normality, homogeneity, linearity
durbinWatsonTest(final); shapiro.test(final$resid)
ncvTest(final); crPlots(lm(PE ~ DE + IDE + VE + CC+ EC + Gender + Age + Residence + Location, data = df))
plot(final)

###### Mixed linear model##########
library(ggplot2)
library(lme4) 
library(car)
library(MuMIn)

df$Gender <- as.factor(df$Gender)
df$Age <- as.factor(df$Age)
df$Residence <- as.factor(df$Residence)
df$Location <- as.factor(df$Location)
df$School <- as.factor(df$School)
summary(df)
dim(df)

## Fitting a mixed-effects model with lmer
mod<- lmer (PE ~ DE + IDE + VE + CC+ EC + Gender + Age + Residence + Location + (1|School), data=df)
par(mfrow=c(1,2)) 
hist(resid(mod))  
str(mod)

# residuals should be normally distributed
hist(ranef(mod)[[1]][,1])
# random variables (intercepts) should be normally distributed
ranef(mod)
summary(mod) 

## Diagnostics (residual properties)
par(mfrow=c(2,2))
resids <- resid(mod, type='pearson')
length(resids)
qqPlot(resid(mod))

## plot the sqrt of the absolute residuals against fitted values
plot(sqrt(abs(resid(mod)))~ fitted(mod))
lines(lowess(sqrt(abs(resid(mod)))~
               fitted(mod)), col='red')

######## Random effects should be normally distributed
res1 <- ranef(mod)$School[,1]
summary(res1) ## mean ~ 0
hist(res1) ## histogram
abline(v=mean(res1), col='red')
par(mfrow=c(1,1)); qqPlot(res1)
qqPlot(ranef(mod)$School$'(Intercept)')

## use the confint function to conduct bootstraps
confint.result <- confint(mod, method='boot', oldNames=F, nsim = 5000) 
confint.result

## get an R2 type statistic (not equivalent to OLS R-square)
r.squaredGLMM(mod)


##### 2.Pro-nature behavior (PE)#######

fit<-lm(PN~ DE+IDE+VE+DE*IDE+DE*VE+IDE*VE++CC+EC+ Gender+Age+Residence+Location+DE*Residence+IDE*Residence+VE*Residence,data = df)
summary(fit)
vif(fit)
plot(fit)
fit_step <- step(fit) 
summary(fit_step)
vif(fit_step)

## Multiple Regression                                               
## Use the full model as a star

fit<-lm(PN~ DE+IDE+VE+DE*IDE+DE*VE+IDE*VE++CC+EC+ Gender+Age+Residence+Location+DE*Residence+IDE*Residence+VE*Residence,data = df);
summary (fit)

## Update model
fit1 <- step (fit) 
summary(fit1)

# Evaluate model
# multicollinearity

vif(lm(PN~DE+IDE+VE+CC+ EC+Gender+Age+Residence+Location , data = df))
final <- lm(PN ~ DE + IDE + VE + CC+ EC + Gender + Age + Residence + Location, data = df)

# Independence, normality, homogeneity, linearity
durbinWatsonTest(final); shapiro.test(final$resid)
ncvTest(final); crPlots(lm(PN ~ DE + IDE + VE + CC+ EC + Gender + Age + Residence + Location, data = df))
plot(final)

###### Mixed linear model##########
## Fitting a mixed-effects model with lmer
mod<- lmer (PN ~ DE + IDE + VE + CC+ EC + Gender + Age + Residence + Location + (1|School), data=df)
par(mfrow=c(1,2)) 
hist(resid(mod))  
str(mod)

# residuals should be normally distributed
hist(ranef(mod)[[1]][,1])
# random variables (intercepts) should be normally distributed
ranef(mod)
summary(mod) 

## Diagnostics (residual properties)
par(mfrow=c(2,2))
resids <- resid(mod, type='pearson')
length(resids)
qqPlot(resid(mod))

## plot the sqrt of the absolute residuals against fitted values
plot(sqrt(abs(resid(mod)))~ fitted(mod))
lines(lowess(sqrt(abs(resid(mod)))~
               fitted(mod)), col='red')

######## Random effects should be normally distributed
res1 <- ranef(mod)$School[,1]
summary(res1) ## mean ~ 0
hist(res1) ## histogram
abline(v=mean(res1), col='red')
par(mfrow=c(1,1)); qqPlot(res1)
qqPlot(ranef(mod)$School$'(Intercept)')

## use the confint function to conduct bootstraps
confint.result <- confint(mod, method='boot', oldNames=F, nsim = 5000) 
confint.result


## get an R2 type statistic (not equivalent to OLS R-square)
r.squaredGLMM(mod)