#instellen working directory
setwd("C:/Users/smcomper/Documents/Documents/Documents/FWO postdoc/BEPAS_ALECS")
dir()

#inlezen data
library(foreign)
mydata <- read.spss("FINALE_DATA (5).sav",header=T, to.data.frame = TRUE)
View(mydata)
attach(mydata)

install.packages("Rcpp")
library(Rcpp)

#installeren en laden package voor multilevel analyses
install.packages("lme4")
library(lme4)

#installeren en laden package om ook p-waarden te krijgen
install.packages("lmerTest")
library(lmerTest)

#normaal verdeelde data => linear mixed model
null_model<-lmer(acc_lightperday ~  1 +
                   (1 | buurtID),
                 data=mydata, REML=FALSE)
summary(null_model)

ICC.Model<-function(Model.Name) {
  tau.null_model<-as.numeric(lapply(summary(Model.Name)$varcor, diag))
  sigma.null_model <- as.numeric(attr(summary(Model.Name)$varcor, "sc")^2)
  ICC.null_model <- tau.null_model/(tau.null_model+sigma.null_model)
  return(ICC.null_model)
}

ICC.Model(null_model)

#modelopbouw
#1. univariaat
#2. alle significante demograf covariaten + ALLE omgevingsfactoren
#3. 1 voor 1 interactie-effecten met leeftijdsgroep toevoegen
#stap 2 en 3 herhalen voor 500m buffer

#leeftijd
model_age<-lmer(acc_lightperday ~  1 + age_group +
                  (1 | buurtID),
                data=mydata, REML=FALSE)
summary(model_age)

#geslacht
model_sex<-lmer(acc_lightperday ~  1 + geslacht +
                  (1 | buurtID),
                data=mydata, REML=FALSE)
summary(model_sex)

#educatie
model_educ<-lmer(acc_lightperday ~  1 + educ +
                   (1 | buurtID),
                 data=mydata, REML=FALSE)
summary(model_educ)

#buurtSES
model_buurteduc<-lmer(acc_lightperday ~  1 + SESbuurt +
                   (1 | buurtID),
                 data=mydata, REML=FALSE)
summary(model_buurteduc)

#validedagen
model_validdays<-lmer(acc_lightperday ~  1 + acc_vlddays +
                        (1 | buurtID),
                      data=mydata, REML=FALSE)
summary(model_validdays)

#weartime
model_weartime<-lmer(acc_lightperday ~  1 + acc_timeperday +
                       (1 | buurtID),
                     data=mydata, REML=FALSE)
summary(model_weartime)

#entropy
model_entropy500m<-lmer(acc_lightperday ~  1 + entropy500m_all_DEF +
                          (1 | buurtID),
                        data=mydata, REML=FALSE)
summary(model_entropy500m)

model_entropy1km<-lmer(acc_lightperday ~  1 + entropy1km_all_DEF +
                         (1 | buurtID),
                       data=mydata, REML=FALSE)
summary(model_entropy1km)

#residentialdensity
model_resden500m<-lmer(acc_lightperday ~  1 + res_density_500m +
                         (1 | buurtID),
                       data=mydata, REML=FALSE)
summary(model_resden500m)

model_resden1km<-lmer(acc_lightperday ~  1 + res_density_1km +
                        (1 | buurtID),
                      data=mydata, REML=FALSE)
summary(model_resden1km)

#parkdensity
model_park500m<-lmer(acc_lightperday ~  1 + park_density_500m +
                       (1 | buurtID),
                     data=mydata, REML=FALSE)
summary(model_park500m)

model_park1km<-lmer(acc_lightperday ~  1 + park_density_1km +
                      (1 | buurtID),
                    data=mydata, REML=FALSE)
summary(model_park1km)

#transportdensity
model_trans500m<-lmer(acc_lightperday ~  1 + transport_density_500m +
                        (1 | buurtID),
                      data=mydata, REML=FALSE)
summary(model_trans500m)

model_trans1km<-lmer(acc_lightperday ~  1 + transport_density_1km +
                       (1 | buurtID),
                     data=mydata, REML=FALSE)
summary(model_trans1km)

#intersectiondensity
model_intersec500m<-lmer(acc_lightperday ~  1 + intersection_density_500m +
                           (1 | buurtID),
                         data=mydata, REML=FALSE)
summary(model_intersec500m)

model_intersec1km<-lmer(acc_lightperday ~  1 + intersection_density_1km +
                          (1 | buurtID),
                        data=mydata, REML=FALSE)
summary(model_intersec1km)


#model met alle relevante covariaten en omgevingsfactoren 500m
model_multi_500m<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                         park_density_500m + res_density_500m + entropy500m_all_DEF + transport_density_500m +intersection_density_500m +
                         (1 | buurtID),
                       data=mydata, REML=FALSE)
summary(model_multi_500m)

#modellen met interactie-effecten 1 km 
model_multi_500m_intersec1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                   park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                   transport_density_500m +intersection_density_500m + age_group*intersection_density_500m +
                                   (1 | buurtID),
                                 data=mydata, REML=FALSE)
summary(model_multi_500m_intersec1)

anova(model_multi_500m,model_multi_500m_intersec1)

model_multi_500m_intersec2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                   park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                   transport_density_500m +intersection_density_500m + relevel(age_group,ref= 'adult')*intersection_density_500m +
                                   (1 | buurtID),
                                 data=mydata, REML=FALSE)
summary(model_multi_500m_intersec2)

model_multi_500m_residential1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                      park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                      transport_density_500m +intersection_density_500m + age_group*res_density_500m +
                                      (1 | buurtID),
                                    data=mydata, REML=FALSE)
summary(model_multi_500m_residential1)

anova(model_multi_500m,model_multi_500m_residential1)

model_multi_500m_residential2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                      park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                      transport_density_500m +intersection_density_500m + relevel(age_group,ref= 'adult')*res_density_500m +
                                      (1 | buurtID),
                                    data=mydata, REML=FALSE)
summary(model_multi_500m_residential2)

model_multi_500m_entropy1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                  park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                  transport_density_500m +intersection_density_500m + age_group*entropy500m_all_DEF +
                                  (1 | buurtID),
                                data=mydata, REML=FALSE)
summary(model_multi_500m_entropy1)

anova(model_multi_500m,model_multi_500m_entropy1)

model_multi_500m_entropy2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                  park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                  transport_density_500m +intersection_density_500m + relevel(age_group,ref= 'adult')*entropy500m_all_DEF +
                                  (1 | buurtID),
                                data=mydata, REML=FALSE)
summary(model_multi_500m_entropy2)

model_multi_500m_park1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                               park_density_500m + res_density_500m + entropy500m_all_DEF + 
                               transport_density_500m +intersection_density_500m + age_group*park_density_500m +
                               (1 | buurtID),
                             data=mydata, REML=FALSE)
summary(model_multi_500m_park1)

anova(model_multi_500m, model_multi_500m_park1)

model_multi_500m_park2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                               park_density_500m + res_density_500m + entropy500m_all_DEF + 
                               transport_density_500m +intersection_density_500m + relevel(age_group,ref= 'adult')*park_density_500m +
                               (1 | buurtID),
                             data=mydata, REML=FALSE)
summary(model_multi_500m_park2)

model_multi_500m_transport1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                    park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                    transport_density_500m +intersection_density_500m + age_group*transport_density_500m +
                                    (1 | buurtID),
                                  data=mydata, REML=FALSE)
summary(model_multi_500m_transport1)

anova(model_multi_500m, model_multi_500m_transport1)

model_multi_500m_transport2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                    park_density_500m + res_density_500m + entropy500m_all_DEF + 
                                    transport_density_500m +intersection_density_500m + relevel(age_group,ref= 'adult')*transport_density_500m +
                                    (1 | buurtID),
                                  data=mydata, REML=FALSE)
summary(model_multi_500m_transport2)

#model met alle relevante covariaten en omgevingsfactoren 1km
model_multi_1km<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                        park_density_1km + res_density_1km + entropy1km_all_DEF + transport_density_1km +intersection_density_1km +
                        (1 | buurtID),
                      data=mydata, REML=FALSE)
summary(model_multi_1km)

#modellen met interactie-effecten 1 km 
model_multi_1km_intersec1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                  park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                  transport_density_1km +intersection_density_1km + age_group*intersection_density_1km +
                                  (1 | buurtID),
                                data=mydata, REML=FALSE)
summary(model_multi_1km_intersec1)

anova(model_multi_1km, model_multi_1km_intersec1)

model_multi_1km_intersec2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                  park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                  transport_density_1km +intersection_density_1km + relevel(age_group,ref= 'adult')*intersection_density_1km +
                                  (1 | buurtID),
                                data=mydata, REML=FALSE)
summary(model_multi_1km_intersec2)

model_multi_1km_residential1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                     park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                     transport_density_1km +intersection_density_1km + age_group*res_density_1km +
                                     (1 | buurtID),
                                   data=mydata, REML=FALSE)
summary(model_multi_1km_residential1)

anova(model_multi_1km,model_multi_1km_residential1)

model_multi_1km_residential2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                     park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                     transport_density_1km +intersection_density_1km + relevel(age_group,ref= 'adult')*res_density_1km +
                                     (1 | buurtID),
                                   data=mydata, REML=FALSE)
summary(model_multi_1km_residential2)

model_multi_1km_entropy1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                 park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                 transport_density_1km +intersection_density_1km + age_group*entropy1km_all_DEF +
                                 (1 | buurtID),
                               data=mydata, REML=FALSE)
summary(model_multi_1km_entropy1)

anova(model_multi_1km, model_multi_1km_entropy1)

model_multi_1km_entropy2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                 park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                 transport_density_1km +intersection_density_1km + relevel(age_group,ref= 'adult')*entropy1km_all_DEF +
                                 (1 | buurtID),
                               data=mydata, REML=FALSE)
summary(model_multi_1km_entropy2)

model_multi_1km_park1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                              park_density_1km + res_density_1km + entropy1km_all_DEF + 
                              transport_density_1km +intersection_density_1km + age_group*park_density_1km +
                              (1 | buurtID),
                            data=mydata, REML=FALSE)
summary(model_multi_1km_park1)

anova(model_multi_1km, model_multi_1km_park1)

model_multi_1km_park2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                              park_density_1km + res_density_1km + entropy1km_all_DEF + 
                              transport_density_1km +intersection_density_1km + relevel(age_group,ref= 'adult')*park_density_1km +
                              (1 | buurtID),
                            data=mydata, REML=FALSE)
summary(model_multi_1km_park2)

model_multi_1km_transport1<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                   park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                   transport_density_1km +intersection_density_1km + age_group*transport_density_1km +
                                   (1 | buurtID),
                                 data=mydata, REML=FALSE)
summary(model_multi_1km_transport1)

anova(model_multi_1km, model_multi_1km_transport1)

model_multi_1km_transport2<-lmer(acc_lightperday ~  1 + age_group + geslacht + educ + acc_timeperday + 
                                   park_density_1km + res_density_1km + entropy1km_all_DEF + 
                                   transport_density_1km +intersection_density_1km + relevel(age_group,ref= 'adult')*transport_density_1km +
                                   (1 | buurtID),
                                 data=mydata, REML=FALSE)
summary(model_multi_1km_transport2)