Wdata=import("C:/Users/Eason Zhang/Dropbox/social impact project/employeedata.sav")%>%as.data.table()
numerical_variable_names <- names(Wdata)[sapply(Wdata, is.numeric)]
Wdata=group_mean_center(Wdata, numerical_variable_names,by="FactoryAssessedID", add.suffix=".GroC")
BMdata=import("C:/Users/Eason Zhang/Dropbox/social impact project/managerdata.sav")%>%as.data.table()
BMdata1=BMdata[year == 2019]
BMdata2=BMdata[year == 2020]
BMdata_mean <- BMdata[, .SD, .SDcols = sapply(BMdata, is.numeric)]
BMdata_mean <- BMdata_mean[, lapply(.SD, mean, na.rm = TRUE), by = .(FactoryAssessedID)]
BAdata=import("C:/Users/Eason Zhang/Dropbox/social impact project/Auditsdata.sav")%>%as.data.table()
BAdata1=BAdata[year == 2019]
BAdata2=BAdata[year == 2020]
setnames(BMdata1, old = names(BMdata1), new = paste0("BM19.", names(BMdata1)))
setnames(BMdata2, old = names(BMdata2), new = paste0("BM20.", names(BMdata2)))
setnames(BMdata_mean, old = names(BMdata_mean), new = paste0("BM.", names(BMdata_mean)))
setnames(BAdata1, old = names(BAdata1), new = paste0("BA19.", names(BAdata)))
setnames(BAdata2, old = names(BAdata2), new = paste0("BA20.", names(BAdata)))
setnames(BMdata_mean, "BM.FactoryAssessedID", "FactoryAssessedID")
setnames(BMdata1, "BM19.FactoryAssessedID", "FactoryAssessedID")
setnames(BMdata2, "BM20.FactoryAssessedID", "FactoryAssessedID")
setnames(BAdata1, "BA19.FactoryAssessedID", "FactoryAssessedID")
setnames(BAdata2, "BA20.FactoryAssessedID", "FactoryAssessedID")
Bdata<- Reduce(function(x, y) merge(x, y, by = "FactoryAssessedID", all = TRUE),list(BAdata1, BAdata2, BMdata_mean, BMdata1, BMdata2))
BWdata <- merge(Wdata, Bdata, by = "FactoryAssessedID", all = FALSE)
BWdata$clus <- BWdata$FactoryAssessedID
Bdata$clus <- Bdata$FactoryAssessedID
export(BMdata1,"managerdata19.sav")
export(BMdata2,"managerdata20.sav")
export(BMdata_mean,"managerdata_mean.sav")
export(BAdata1,"Auditsdata19.sav")
export(BAdata2,"Auditsdata20.sav")
export(Bdata,"Bdata.sav")
export(BWdata,"BWdata.sav")
Illu.data=data
# Assuming Illu.data is already a data.table; if not, convert it using setDT(Illu.data)
setnames(Illu.data, old = c("healthown", "migrant10.GroC", "migrant10", "migrant10_mean", "BA19.complall", "BM.responsiblesourcing"),
new = c("Y","Xw", "X", "Xb", "M", "W"))
Illu.Bdata=Bdata
# Assuming Illu.data is already a data.table; if not, convert it using setDT(Illu.data)
setnames(Illu.Bdata, old = c("BA19.complall", "BM.responsiblesourcing"),
new = c("M", "W"))
MLM.ab2=lmer(M~W + (1|clus), na.action = na.exclude, data = Illu.data, control=lmerControl(optimizer="bobyqa")) # Xb is not included for comparing with OLR.ab2
#MLM.ab2=lmer(BA19.complall~BM.responsiblesourcing+migrant10_mean + (1|clus), na.action = na.exclude, data = data, control=lmerControl(optimizer="bobyqa"))
#HLM_summary(MLM.ab2)
MLM.cw2=lmer(Y~Xb*M + Xb*W + X*W + X*M + (X|clus), na.action = na.exclude, data = data, control=lmerControl(optimizer="bobyqa"))
#MLM.cw2=lmer(Y~X+migrant10_mean + BA19.complall + BM.responsiblesourcing + migrant10_mean:BA19.complall + migrant10_mean:BM.responsiblesourcing + migrant10*BM.responsiblesourcing + migrant10*BA19.complall+ (migrant10|clus), na.action = na.exclude, data = data, control=lmerControl(optimizer="bobyqa"))
#Full version according to statistical model: X + Xb + M + W + Xb:M + Xb:W + X*W + X*M + (X|clus)
#HLM_summary(MLM.cw2)
interact_plot(MLM.cw2, pred = X, modx = M,#Basic setup
modx.values = "plus-minus", modx.labels= c("Low group", "High group"),legend.main="BA19.complall",)+#Set moderators in plot
ylab("Health")+xlab("International immigrant")#+#Set labels of X and Y
## [1] "(Intercept)" "Xb" "M" "W" "X"
## [6] "Xb:M" "Xb:W" "W:X" "M:X"
##
## Model Summary
##
## ────────────────────────────────────────────────────
## (1) Y (2) M
## ────────────────────────────────────────────────────
## (Intercept) 4.739 ** 0.880 ***
## (1.495) (0.002)
## Xb -1.316
## (2.292)
## M -0.909
## (1.680)
## W 0.078 0.003 ***
## (0.061) (0.001)
## X -0.815
## (1.103)
## Xb:M 1.682
## (2.588)
## Xb:W -0.025
## (0.086)
## W:X -0.059
## (0.050)
## M:X 1.121
## (1.255)
## ────────────────────────────────────────────────────
## Marginal R^2 0.012 0.330
## Conditional R^2 0.066 1.000
## AIC 8611.505 -87596.577
## BIC 8690.787 -87572.183
## Num. obs. 3290 3290
## Num. groups: clus 51 51
## Var: clus (Intercept) 0.087 0.000
## Var: clus X 0.141
## Cov: clus (Intercept) X -0.102
## Var: Residual 0.772 0.000
## ────────────────────────────────────────────────────
## Note. * p < .05, ** p < .01, *** p < .001.
##
## Effects of Mediation (rep=5000, seed=1223 if not preset):
## ────────────────────────────────────────────────────────────────────────────────────
## Effects ab SEab a SEa b SEb abLLCI abULCI sig
## ────────────────────────────────────────────────────────────────────────────────────
## 1 Monte Carlo based Mediation 0.006 0.009 1.682 2.588 0.003 0.001 -0.011 0.024 no
## 2 Monte Carlo based Mediation 0.004 0.004 1.121 1.255 0.003 0.001 -0.004 0.013 no
## ────────────────────────────────────────────────────────────────────────────────────
OLS.ab2 <- lm(M ~ W, data = Illu.data, na.action = na.exclude)
OLS.Me=AutoMC.Me(MLM.cw2,OLS.ab2,cc("Xb:M,M:X"))#,note=T,Cnote=T)#,cc("M"))
##
## Model Summary
##
## ──────────────────────────────────────────────────
## (1) Y (2) M
## ──────────────────────────────────────────────────
## (Intercept) 4.739 ** 0.878 ***
## (1.495) (0.002)
## Xb -1.316
## (2.292)
## M -0.909
## (1.680)
## W 0.078 0.008 ***
## (0.061) (0.001)
## X -0.815
## (1.103)
## Xb:M 1.682
## (2.588)
## Xb:W -0.025
## (0.086)
## W:X -0.059
## (0.050)
## M:X 1.121
## (1.255)
## ──────────────────────────────────────────────────
## Marginal R^2 0.012
## Conditional R^2 0.066
## AIC 8611.505
## BIC 8690.787
## Num. obs. 3290 3290
## Num. groups: clus 51
## Var: clus (Intercept) 0.087
## Var: clus X 0.141
## Cov: clus (Intercept) X -0.102
## Var: Residual 0.772
## R^2 0.064
## Adj. R^2 0.064
## ──────────────────────────────────────────────────
## Note. * p < .05, ** p < .01, *** p < .001.
##
## Effects of Mediation (rep=5000, seed=1223 if not preset):
## ────────────────────────────────────────────────────────────────────────────────────
## Effects ab SEab a SEa b SEb abLLCI abULCI sig
## ────────────────────────────────────────────────────────────────────────────────────
## 1 Monte Carlo based Mediation 0.014 0.021 1.682 2.588 0.008 0.001 -0.028 0.056 no
## 2 Monte Carlo based Mediation 0.009 0.010 1.121 1.255 0.008 0.001 -0.011 0.030 no
## ────────────────────────────────────────────────────────────────────────────────────
combinedResults <- Combine.AutoMC(
autoMCList = list(OLS.Me, MLM.Me),
titles = c("MEDIATION RESULTS BASED ON OLS", "MEDIATION RESULTS BASED ON MLM"))
## ───────────────────────────────────────────────────────────────────────────────────────
## Effects ab SEab a SEa b SEb abLLCI abULCI sig
## ───────────────────────────────────────────────────────────────────────────────────────
## 1 MEDIATION RESULTS BASED ON OLS
## 2 Monte Carlo based Mediation 0.014 0.021 1.682 2.588 0.008 0.001 -0.028 0.056 no
## 3 Monte Carlo based Mediation 0.009 0.010 1.121 1.255 0.008 0.001 -0.011 0.030 no
## 4 MEDIATION RESULTS BASED ON MLM
## 5 Monte Carlo based Mediation 0.006 0.009 1.682 2.588 0.003 0.001 -0.011 0.024 no
## 6 Monte Carlo based Mediation 0.004 0.004 1.121 1.255 0.003 0.001 -0.004 0.013 no
## ───────────────────────────────────────────────────────────────────────────────────────
OLS.ce=ESCaculatorMC(params = OLS.Me$ab[1:2], SEs = OLS.Me$SEab[1:2], expression = "p1-p2")#,note=T,Cnote=T)
## Parameters and SE:
## ──────────────────────
## p1 p2
## ──────────────────────
## Value 0.01375 0.00919
## SE 0.02113 0.01026
## ──────────────────────
## Monte Carlo based Estimation (rep=5000, seed=1223 if not preset):
## ────────────────────────────────────────────────────────────────────
## Expression Estimate MCEstimate S.E. LLCI ULCI sig
## ────────────────────────────────────────────────────────────────────
## Monte Carlo p1-p2 0.005 0.005 (0.023) -0.041 0.051 no
## ────────────────────────────────────────────────────────────────────
## Parameters and SE:
## ────────────────────────
## p1 p2
## ────────────────────────
## Value 0.005578 0.003732
## SE 0.008780 0.004287
## ────────────────────────
## Monte Carlo based Estimation (rep=5000, seed=1223 if not preset):
## ────────────────────────────────────────────────────────────────────
## Expression Estimate MCEstimate S.E. LLCI ULCI sig
## ────────────────────────────────────────────────────────────────────
## Monte Carlo p1-p2 0.002 0.002 (0.010) -0.017 0.021 no
## ────────────────────────────────────────────────────────────────────
Full results: https://www.dropbox.com/scl/fi/u3apr3ncp15qwy6qvfr9h/JAPanalysis1.xlsx?rlkey=c3rhfpyp67jt80b2o7vhsg1lw&dl=0
##-Running 1-1000
###-A. Writing formula
source("https://www.dropbox.com/scl/fi/cgl95x7kbel2r4nqbg9wh/MAIN-EFFECT-ab1-bb2.R?rlkey=39mf8r87r5uu8e3jqgpz5oy2n&raw=1")
###B-E. Run formulas
source("https://www.dropbox.com/scl/fi/628ipoqcn2dia89w74x45/lmerRun.EN.r?rlkey=duk4gv7dognek4p1owuhi1eq9&raw=1")
###-F. Save excel
write.xlsx(Output.lmer, file="JAPanalysis1.xlsx", sheetName="Main.ab1bb2", append=TRUE, row.names=TRUE)
##-Running 1-1000
###-A. Writing formula
source("https://www.dropbox.com/scl/fi/fb7bcj5bwlnuuzi6x80jm/olsMAIN-EFFECT-ab1-bb2.R?rlkey=zvvlr9s18xw66k8vtgzooi040&raw=1")
###B-E. Run formulas
source("https://www.dropbox.com/scl/fi/cbnzb2oyqy5667s368cgt/lmRun.r?rlkey=wsem22it5sj6ta3ii3vne237i&raw=1")
###-F. Save excel
write.xlsx(Output.lmer, file="JAPanalysis1.xlsx", sheetName="olsMain.ab1bb2", append=TRUE, row.names=TRUE)