##################################################################################
##### GLMM analyses of the PMs effects on AD symptoms by SSC weather type ########
##################################################################################
library(lme4)

SEL <- read.csv(file="C://Users//Eung Ki Jung//Desktop//SSC and AD//AD_PMs_twa_seoul_neat.csv", as.is=T)

ssc.1 <- subset(SEL, ssc =="1")
ssc.2 <- subset(SEL, ssc =="2")
ssc.3 <- subset(SEL, ssc =="3")
ssc.4 <- subset(SEL, ssc =="4")
ssc.5 <- subset(SEL, ssc =="5")
ssc.6 <- subset(SEL, ssc =="6")
ssc.66 <- subset(SEL, ssc =="66")
ssc.7 <- subset(SEL, ssc =="7")
ssc.6 <- rbind(ssc.6, ssc.66)


PM25 <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
PM25.dm <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
PM25.dp <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.2, binomial, verbose =2L)
PM25.dt <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.3, binomial, verbose =2L)
PM25.mm <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.4, binomial, verbose =2L)
PM25.mp <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.5, binomial, verbose =2L)
PM25.mt <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.6, binomial, verbose =2L)
PM25.tr <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.7, binomial, verbose =2L)

coef.pm25=summary(PM25)$coefficients[2,]
coef.pm25.dm=summary(PM25.dm)$coefficients[2,]
coef.pm25.dp=summary(PM25.dp)$coefficients[2,]
coef.pm25.dt=summary(PM25.dt)$coefficients[2,]
coef.pm25.mm=summary(PM25.mm)$coefficients[2,]
coef.pm25.mp=summary(PM25.mp)$coefficients[2,]
coef.pm25.mt=summary(PM25.mt)$coefficients[2,] 
coef.pm25.tr=summary(PM25.tr)$coefficients[2,]
coef <- rbind(coef.pm25, coef.pm25.dm, coef.pm25.dp,coef.pm25.dt,coef.pm25.mm, coef.pm25.mp, coef.pm25.mt, coef.pm25.tr)
coef

OR=exp(12.72*coef[,1])
SE=12.72*coef[,2]
Pvalue=coef[,4]
ORhi=exp(12.72*coef[,1]+12.72*1.96*coef[,2])
ORlo=exp(12.72*coef[,1]-12.72*1.96*coef[,2])
OR.rng <- rbind(OR,SE,ORlo,ORhi, Pvalue)
colnames(OR.rng) <-c("PM2.5","DM","DP","DT","MM","MP","MT","Tr")
t(OR.rng)

fit0 <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(q1_7) + factor(q4) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit01 <- glmer(Atopy ~ PM25_ta_lag1 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit02 <- glmer(Atopy ~ PM25_ta_lag2 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit03 <- glmer(Atopy ~ PM25_ta_lag3 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit04 <- glmer(Atopy ~ PM25_ta_lag4 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit05 <- glmer(Atopy ~ PM25_ta_lag5 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)

fit1 <- glmer(Atopy ~ PM25_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(q1_7) + factor(q4) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit11 <- glmer(Atopy ~ PM25_ta_lag1 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit12 <- glmer(Atopy ~ PM25_ta_lag2 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit13 <- glmer(Atopy ~ PM25_ta_lag3 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit14 <- glmer(Atopy ~ PM25_ta_lag4 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit15 <- glmer(Atopy ~ PM25_ta_lag5 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)

PM25_lag0 =summary(fit0)$coefficients[2,]
PM25_lag1 =summary(fit01)$coefficients[2,]
PM25_lag2 =summary(fit02)$coefficients[2,]
PM25_lag3 =summary(fit03)$coefficients[2,]
PM25_lag4 =summary(fit04)$coefficients[2,]
PM25_lag5 =summary(fit05)$coefficients[2,]
coef <- rbind(PM25_lag0, PM25_lag1,PM25_lag2,PM25_lag3, PM25_lag4, PM25_lag5)
coef

OR=exp(12.72*coef[,1])
SE=12.72*coef[,2]
Pvalue=coef[,4]
ORhi=exp(12.72*coef[,1]+12.72*1.96*coef[,2])
ORlo=exp(12.72*coef[,1]-12.72*1.96*coef[,2])
OR.rng <- rbind(OR,SE,ORlo,ORhi, Pvalue)
colnames(OR.rng) <-c("lag0","lag1","lag2","lag3","lag4","lag5")
t(OR.rng)

DM_lag0 =summary(fit1)$coefficients[2,]
DM_lag1 =summary(fit11)$coefficients[2,]
DM_lag2 =summary(fit12)$coefficients[2,]
DM_lag3 =summary(fit13)$coefficients[2,]
DM_lag4 =summary(fit14)$coefficients[2,]
DM_lag5 =summary(fit15)$coefficients[2,]
coef <- rbind(DM_lag0, DM_lag1,DM_lag2, DM_lag3, DM_lag4, DM_lag5)
coef

OR=exp(12.72*coef[,1])
SE=12.72*coef[,2]
Pvalue=coef[,4]
ORhi=exp(12.72*coef[,1]+12.72*1.96*coef[,2])
ORlo=exp(12.72*coef[,1]-12.72*1.96*coef[,2])
OR.rng <- rbind(OR,SE,ORlo,ORhi, Pvalue)
colnames(OR.rng) <-c("lag0","lag1","lag2","lag3","lag4","lag5")
t(OR.rng)

##### PM10_TWA analyses ##########
PM10 <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(q1_7) + factor(q4) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
PM10.dm <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
PM10.dp <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.2, binomial, verbose =2L)
PM10.dt <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.3, binomial, verbose =2L)
PM10.mm <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.4, binomial, verbose =2L)
PM10.mp <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.5, binomial, verbose =2L)
PM10.mt <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.6, binomial, verbose =2L)
PM10.tr <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.7, binomial, verbose =2L)

coef.pm10=summary(PM10)$coefficients[2,]
coef.pm10.dm=summary(PM10.dm)$coefficients[2,]
coef.pm10.dp=summary(PM10.dp)$coefficients[2,]
coef.pm10.dt=summary(PM10.dt)$coefficients[2,]
coef.pm10.mm=summary(PM10.mm)$coefficients[2,]
coef.pm10.mp=summary(PM10.mp)$coefficients[2,]
coef.pm10.mt=summary(PM10.mt)$coefficients[2,] 
coef.pm10.tr=summary(PM10.tr)$coefficients[2,]
coef <- rbind(coef.pm10, coef.pm10.dm, coef.pm10.dp,coef.pm10.dt,coef.pm10.mm, coef.pm10.mp, coef.pm10.mt, coef.pm10.tr)
coef
OR=exp(23.99*coef[,1])
SE=23.99*coef[,2]
Pvalue=coef[,4]
ORhi=exp(23.99*coef[,1]+23.99*1.96*coef[,2])
ORlo=exp(23.99*coef[,1]-23.99*1.96*coef[,2])
OR.rng <- rbind(OR,ORlo,ORhi, Pvalue)
colnames(OR.rng) <-c("PM10","DM","DP","DT","MM","MP","MT","Tr")
t(OR.rng)

fit2 <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow) + factor(season)+ SCORAD  + factor(q1_7) + factor(q4) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit21 <- glmer(Atopy ~ PM10_ta_lag1 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit22 <- glmer(Atopy ~ PM10_ta_lag2 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit23 <- glmer(Atopy ~ PM10_ta_lag3 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit24 <- glmer(Atopy ~ PM10_ta_lag4 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)
fit25 <- glmer(Atopy ~ PM10_ta_lag5 + Tem + RH + age + factor(sex)+ factor(dow)+ factor(season) + SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=SEL, binomial, verbose =2L)

fit3 <- glmer(Atopy ~ PM10_ta + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(q1_7) + factor(q4) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit31 <- glmer(Atopy ~ PM10_ta_lag1 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit32 <- glmer(Atopy ~ PM10_ta_lag2 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit33 <- glmer(Atopy ~ PM10_ta_lag3 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit34 <- glmer(Atopy ~ PM10_ta_lag4 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)
fit35 <- glmer(Atopy ~ PM10_ta_lag5 + Tem + RH + age + factor(sex)+ factor(dow)+ SCORAD  + factor(fever) + factor(TCS) + (1|ANGEL_ID), data=ssc.1, binomial, verbose =2L)

PM10_lag0 =summary(fit2)$coefficients[2,]
PM10_lag1 =summary(fit21)$coefficients[2,]
PM10_lag2 =summary(fit22)$coefficients[2,]
PM10_lag3 =summary(fit23)$coefficients[2,]
PM10_lag4 =summary(fit24)$coefficients[2,]
PM10_lag5 =summary(fit25)$coefficients[2,]
coef <- rbind(PM10_lag0,PM10_lag1, PM10_lag2, PM10_lag3, PM10_lag4, PM10_lag5)
coef

OR=exp(23.99*coef[,1])
SE=23.99*coef[,2]
Pvalue=coef[,4]
ORhi=exp(23.99*coef[,1]+23.99*1.96*coef[,2])
ORlo=exp(23.99*coef[,1]-23.99*1.96*coef[,2])
OR.rng <- rbind(OR,ORlo,ORhi, Pvalue)
colnames(OR.rng) <-c("lag0", "lag1","lag2","lag3","lag4","lag5")
t(OR.rng)

PM10_lag0_dm =summary(fit3)$coefficients[2,]
PM10_lag1_dm =summary(fit31)$coefficients[2,]
PM10_lag2_dm =summary(fit32)$coefficients[2,]
PM10_lag3_dm =summary(fit33)$coefficients[2,]
PM10_lag4_dm =summary(fit34)$coefficients[2,]
PM10_lag5_dm =summary(fit35)$coefficients[2,]
coef <- rbind(PM10_lag0_dm,PM10_lag1_dm, PM10_lag2_dm, PM10_lag3_dm, PM10_lag4_dm, PM10_lag5_dm)
coef


OR=exp(23.99*coef[,1])
SE=23.99*coef[,2]
Pvalue=coef[,4]
ORhi=exp(23.99*coef[,1]+23.99*1.96*coef[,2])
ORlo=exp(23.99*coef[,1]-23.99*1.96*coef[,2])
OR.rng <- rbind(OR,ORlo,ORhi, Pvalue)
colnames(OR.rng) <-c("lag0", "lag1","lag2","lag3","lag4","lag5")
t(OR.rng)