library(haven)
datasekolah <- read_sav(file = "C:/datapm/pisa2012_sch_IDN.sav")
datasiswa <- read_sav(file = "C:/datapm/pisa2012_stud_IDN.sav") 

#menggabungkan variabel perhatian
d1 <- data.frame(datasiswa[,c("SCHOOLID","ST04Q01","ST73Q01","ST73Q02","ST74Q01","ST74Q02",
                              "ST75Q01","ST75Q02","ST76Q01","ST76Q02","PV1MAPE","PV2MAPE",
                              "PV3MAPE","PV4MAPE","PV5MAPE","PV1MAPF","PV2MAPF","PV3MAPF",
                              "PV4MAPF","PV5MAPF","PV1MAPI","PV2MAPI","PV3MAPI","PV4MAPI",
                              "PV5MAPI","CLSMAN")])

#Data Level 2 
d2 <- data.frame(datasekolah[,c("SCHOOLID","SC13Q01","SC13Q02","SC13Q03")])

#menggabungkan variabel perhatian
library(dplyr)
d3 <- full_join(d1,d2, by=c("SCHOOLID"))

#hapus baris dengan NA
d4 <- d3[complete.cases(d3),]

#================================================================================================================
d4%>%count(SCHOOLID)

which(d4$SCHOOLID=='0000001')
which(d4$SCHOOLID=='0000009')
which(d4$SCHOOLID=='0000035')
which(d4$SCHOOLID=='0000054')
which(d4$SCHOOLID=='0000104')
which(d4$SCHOOLID=='0000135')
which(d4$SCHOOLID=='0000143')
which(d4$SCHOOLID=='0000146')

d5 <- d4[-c(1:2,47,206:207,308:309,670:671,915:916,979:980,986:987),]

#=================================================================================================================

#Gender: 1 Female 2 Male => 1:Male 0:Female
d5$ST04Q01<- factor(d5$ST04Q01)
d5$Gender <- factor(ifelse(d5$ST04Q01=="2",1,0))

#Internet School Work
d5$ISW <- rowSums(d5[,c("SC13Q01","SC13Q02","SC13Q03")])

#OTL
d5$OTL1 <- ifelse(d5$ST73Q01 == 4, 1,
                  + ifelse(d5$ST73Q01 == 3, 2,
                           + ifelse(d5$ST73Q01 == 2, 3, 4)))
d5$OTL2 <- ifelse(d5$ST73Q02 == 4, 1,
                 + ifelse(d5$ST73Q02 == 3, 2,
                          + ifelse(d5$ST73Q02 == 2, 3, 4)))
d5$OTL3 <- ifelse(d5$ST74Q01 == 4, 1,
                   + ifelse(d5$ST74Q01 == 3, 2,
                            + ifelse(d5$ST74Q01 == 2, 3, 4)))
d5$OTL4 <- ifelse(d5$ST74Q02 == 4, 1,
                  + ifelse(d5$ST74Q02 == 3, 2,
                           + ifelse(d5$ST74Q02 == 2, 3, 4)))
d5$OTL5 <- ifelse(d5$ST75Q01 == 4, 1,
                   + ifelse(d5$ST75Q01 == 3, 2,
                            + ifelse(d5$ST75Q01 == 2, 3, 4)))
d5$OTL6 <- ifelse(d5$ST75Q02 == 4, 1,
                  + ifelse(d5$ST75Q02 == 3, 2,
                           + ifelse(d5$ST75Q02 == 2, 3, 4)))
d5$OTL7 <- ifelse(d5$ST76Q01 == 4, 1,
                   + ifelse(d5$ST76Q01 == 3, 2,
                            + ifelse(d5$ST76Q01 == 2, 3, 4)))
d5$OTL8 <- ifelse(d5$ST76Q02 == 4, 1,
                  + ifelse(d5$ST76Q02 == 3, 2,
                           + ifelse(d5$ST76Q02 == 2, 3, 4)))
d5$OTL <- rowSums(d5[,c("OTL1","OTL2","OTL3","OTL4","OTL5","OTL6","OTL7","OTL8")])

#Nilai Literasi Matematis
d5$EMPLOY <- rowMeans(d5[,c("PV1MAPE","PV2MAPE","PV3MAPE","PV4MAPE","PV5MAPE")])
d5$FORMULATE <- rowMeans(d5[,c("PV1MAPF","PV2MAPF","PV3MAPF","PV4MAPF","PV5MAPF")])
d5$INTERPRET <- rowMeans(d5[,c("PV1MAPI","PV2MAPI","PV3MAPI","PV4MAPI","PV5MAPI")])
d5$mathliteracy <- rowSums(d5[,c("EMPLOY","FORMULATE","INTERPRET")])

#===========================================================================================

library(psych)
describe(d5)

#===========================================================================================

#Model Null
library(nlme)
Model0 <- lme(fixed=mathliteracy~1, random=~1|SCHOOLID, data=d5)
summary(Model0)


#ICC
ICC=(144.5797^2/(144.5797^2+131.9705^2))
ICC

#Design Effect
nc=1473/163
nc
DE=1+(nc-1)*ICC
DE

#Cara lain
library(lme4)
Model0.1 <- lmer(mathliteracy~1+(1|SCHOOLID), REML=FALSE, data=d4)
library(jtools)
summ(Model0.1)

#===================================================
#Group Mean Centering
d5$gc_OTL <- d5$OTL-mean(d5$OTL,na.rm=TRUE)

#Grand Mean Centering
d5$gd_ISW <- d5$ISW-mean(d5$ISW,na.rm=TRUE)

#====================================================
 
#Level 1 fixed
library(nlme)
Model1 <- lme (fixed = mathliteracy~Gender+gc_OTL, random=~1|SCHOOLID, data = d5)
summary(Model1)

#Model level 1 dan 2 fixed
library(nlme)
Model2 <- lme (fixed = mathliteracy~Gender+gc_OTL+gd_ISW, random=~1|SCHOOLID, data = d5)
summary(Model2)

#Model Interaksi
library(nlme)
Model3 <- lme (fixed = mathliteracy~Gender+gc_OTL+gd_ISW+gc_OTL*gd_ISW, random=~gc_OTL|SCHOOLID, data = d5)
summary(Model3)




#Model level 1 dan 2 random
library(nlme)
Model1.2 <- lme (fixed = mathliteracy~gd_OTL+Gender+gd_ISW, random=~gd_OTL|SCHOOLID, data = d5)
summary(Model1.2)

library(lme4)
Model2 <- lmer(mathliteracy~Gender+OTL+ISW+(1|SCHOOLID), REML=FALSE, data=d4)
library(jtools)
summ(Model2)


library(psych)
describe(d4)