#uploading library
library(tidyverse)
library(data.table)
library(readxl)
library(tibble)
library(purrr)
library(mlogit)
library(car)
library(ggplot2)
library(nnet)
library(MASS)
library(Rmisc)
library(sjPlot)
#for multilevel analysis
library(car)
library(ggplot2)
library(lme4)
library(reshape)
library(ez)
library(nlme)
library(rcompanion)


#no outliers - preparing dataset
setwd("C:/Users/Utente/Documents/Matteo PhD/study3/creating dataset PA study 3")

dat <- read.csv("study3datasetPAFINALnooutlier.csv")

dat$classid
dat$sex <- as.factor(dat$Sex)
dat$group <- as.factor(dat$Group)
dat$Ethnicity.code <- as.factor(dat$Ethnicity.code)
dat$Deprivation.Decile<- as.factor(dat$Deprivation.Decile)
dat$Ethnicity.code
table(dat$Ethnicity.code)
dat$teacher <- as.factor(dat$teacher)
dat$time <- dat$time -1
dat<-dat%>%
  mutate(content = replace(x, x %in% c("101", "102", "105","106","110","111","112","113","114","115"), "Dance")) %>%
  mutate(content = replace(content, content %in% c("103", "108", "109", "207","3034","307","308","309"), "aLocomotor")) %>%
  mutate(content = replace(content, content %in% c("107", "203","204", "208", "209", "301","302","305","306","310","311","312","313","314","315"), "BallGames")) %>%
  mutate(content = replace(content, content %in% c("201", "202", "205","206","210","211","212","213","214","215"), "Gymnastic"))%>%
  mutate(outdoors = replace(x, x %in% c("101", "102", "103","105","106","107","108","109","110","111","112","113","114","115","201", "202", "203","204","205","206","207","208","209","210","211","212","213","214","215","308","309"), "0")) %>%
  mutate(outdoors = replace(outdoors, outdoors %in% c("301", "302", "3034","305","306","307","310","311","312","313","314","315"), "1"))
dat$content <- as.factor(dat$content)
dat$outdoors <- as.factor(dat$outdoors)
dat$id <- as.factor(dat$id)


#creating data in percentages and minutes

dat$MVPAperc<- dat$MVPA/dat$lessdur*100
dat$MVPA
dat$MVPAperc
max(dat$MVPAperc)
min(dat$MVPAperc)
max(dat$MVPA)
min(dat$MVPA)
dat$MVPA<- dat$MVPA/60
dat$lessdur<- dat$lessdur/60


#models for paper
m1 <- lmer(MVPA~ group+
             (1|id) ,
           data=dat,REML=TRUE, na.action = na.omit)

m2 <- lmer(MVPA~ group+ sex + Decimal.Age + outdoors + content + lessdur+
             (1|schoolid/id) +  (1|teacher),
           data=dat, REML=TRUE, na.action = na.omit)


# NONLINEAR as reference models for paper
#Compare nonlinear to others
dat$group <- as.numeric(dat$group)
dat<-dat%>%
  mutate(groupNONLINEAR = replace(group, group %in% c(1), 4))
dat$groupNONLINEAR <- as.factor(dat$groupNONLINEAR)
dat$groupNONLINEAR
m3 <- lmer(MVPA~ groupNONLINEAR+
           (1|id) ,
           data=dat,REML=TRUE, na.action = na.omit)

m4 <- lmer(MVPA~ groupNONLINEAR+ sex + Decimal.Age+ outdoors + content + lessdur+
             (1|schoolid/id) +  (1|teacher),
           data=dat,REML=TRUE, na.action = na.omit)


#multilevel models output
tab_model(m1,m2)
tab_model(m3,m4)