Wdata=import("C:/Users/Eason Zhang/Dropbox/social impact project/data/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")
#- Within variables
Wdata=import("C:/Users/Eason Zhang/Dropbox/social impact project/data/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")
#-- Within data check
matching_variables.W <- cc("mentalhealth,stress,wage,satoverall,migrant10.GroC")
contains_variables.W <- names(Wdata) %in% matching_variables.W
matching_variables.W <- names(Wdata)[contains_variables.W]
matching_variables.W
#- Between variables
BA2r.1521=import("C:/Users/Eason Zhang/Dropbox/social impact project/data/BAtargetd vs uni practices 2019 to 2021.dta")%>%as.data.table()
Bdata<- import("Bdata.sav")%>%as.data.table()
Bdata2 <- merge(Bdata, BA2r.1521, by = "FactoryAssessedID", all = T)
#-- Between data check
names(BA2r.1521)
names(Bdata)
T_variables.B <- cc("BM.responsiblesourcing,BM19.responsiblesourcing,BM20.responsiblesourcing,BM.adversesourcing,BM19.adversesourcing,BM20.adversesourcing,BA19discriminateall,BA19unipractices,BA19unipracticeAL,BA19unipracticeslong,BA19hrpractices,BA19unipay,BA19voicepractices,BA19firepractices,BA19worktimepractices,BA19targetedpractices,BA19discriminateall,BA19discriminate4items,BA20unipractices,BA20unipracticeAL,BA20unipracticeslong,BA20hrpractices,BA20unipay,BA20voicepractices,BA20firepractices,BA20worktimepractices,BA20targetedpractices,BA20discriminateall,BA20discriminate4items,BA1920unipractices,BA1920unipracticeAL,BA1920unipracticeslong,BA1920hrpractices,BA1920unipay,BA1920voicepractices,BA1920firepractices,BA1920worktimepractices,BA1920targetedpractices,BA1920discriminateall,BA1920discriminate4items,BA2021unipractices,BA2021unipracticeAL,BA2021unipracticeslong,BA2021hrpractices,BA2021unipay,BA2021voicepractices,BA2021firepractices,BA2021worktimepractices,BA2021targetedpractices,BA2021discriminateall,BA2021discriminate4items,BA19migrantsupc,BA20migrantsupc,BA1920migrantsupc,BA2021migrantsupc")
contains_variables.B <- names(Bdata2) %in% T_variables.B
matching_variables.B <- names(Bdata2)[contains_variables.B]
matching_variables.B
matching_B <- data.table(Mvariable = matching_variables.B)
matching_B$Mv=matching_B$Mvariable
T_B <- data.table(Tvariable = T_variables.B)
merged_data <- merge(matching_B, T_B, by.x = "Mvariable", by.y = "Tvariable", all = T)#%>%print_table()
merged_data <- merged_data[order(Mv)]%>%print_table()
#BA targeted uni practices 15 21.dta没有unipracticesAL这个变量
BA2r.1521=import("C:/Users/Eason Zhang/Dropbox/social impact project/data/BA targeted uni practices 15 21.dta")%>%as.data.table()
targeted_variables <- c("unipractices", "unipracticesAL", "unipracticeslong",
"hrpractices", "unipay", "voicepractices",
"firepractices", "worktimepractices", "targetedpractices",
"discriminateall")
contains_variables <- names(BA2r.1521) %in% targeted_variables
matching_variables.1521 <- names(BA2r.1521)[contains_variables]
matching_df.1521 <- data.table(variable = matching_variables.1521)
targeted_df.1521 <- data.table(variable = targeted_variables)
matching_dt.1521 <- matching_df.1521[order(variable)]
targeted_dt.1521 <- targeted_df.1521[order(variable)]
print_table(matching_dt.1521)
print_table(targeted_dt.1521)
#- Full data
BWdata2 <- merge(Bdata2, Wdata, by = "FactoryAssessedID", all = T)
#- Export
export(Bdata2,"Bdata2.sav")
export(BWdata2,"BWdata2.sav")
There is a lot of missing data in our database.
dat=import("BWdata2.sav")%>%as.data.table()
dat$clus <- dat$FactoryAssessedID
#data=rio::import("https://www.dropbox.com/scl/fi/1vhn15jla3duzb14z3r7y/BWdata.sav?rlkey=ek6d7kfw4mopnc7mc20g6nndn&raw=1")%>%as.data.table()
#names(data)
#Freq(data$migrant10.GroC)
#Freq(dat$round)
Bdata<- import("Bdata2.sav")%>%as.data.table()
#names(Bdata)
#Freq(Bdata$FactoryAssessedID)
Freq(Bdata$BM20.B3_1)
## Frequency Statistics:
## ──────────────
## N %
## ──────────────
## (NA) 92 100.0
## ──────────────
## Total N = 92
## Valid N = 0
#HLM_ICC_rWG(data, group="clus", icc.var="newskillpc")
# Assuming 'data' is your data frame
# Find variable names that contain 'childorforcedlabor'
matching_vars <- grep("FOA", names(data), value = TRUE)
# Print the matching variable names
print(matching_vars)
## character(0)
Illu.data=dat[round == 3]
# Assuming Illu.data is already a data.table; if not, convert it using setDT(Illu.data)
setnames(Illu.data, old = c("stress", "migrant10.GroC", "migrant10", "migrant10_mean", "BA20.discrimination", "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(BA19unipractices~BM.responsiblesourcing+ (1|clus), na.action = na.exclude, data = Illu.data, control=lmerControl(optimizer="bobyqa")) # Xb is not
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)
#COEFdata=coef(summary(MLM.ab2))
#View(COEFdata)
#MLM.ab2=lmer(BA19.localworkerpc~BA19.complnoFOA+ (1|clus), na.action = na.exclude, data = data, control=lmerControl(optimizer="bobyqa"))
#HLM_summary(MLM.ab2)
#model_summary(MLM.ab2)
#MLM.cw2=lmer(Y~Xb + X*W + X*M + (X|clus), na.action = na.exclude, data = Illu.data, control=lmerControl(optimizer="bobyqa"))
#MLM.cw2=lmer(Y ~ X*M + (X|clus), na.action = na.exclude, data = Illu.data, control=lmerControl(optimizer="bobyqa"))
MLM.cw2=lmer(Y~Xb*M + Xb*W + X*W + X*M + (X|clus), na.action = na.exclude, data = Illu.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="BA20.discrimination",)+#Set moderators in plot
ylab("stress")+xlab("International immigrant")#+#Set labels of X and Y
## [1] "(Intercept)" "Xb" "M" "W" "X"
## [6] "Xb:M" "Xb:W" "W:X" "M:X"
#MLM.Me=AutoMC.Me(MLM.cw2,MLM.ab2,cc("X:M"))#,note=T,Cnote=T)#,cc("M"))
MLM.Me=AutoMC.Me(MLM.cw2,MLM.ab2,cc("Xb:M,M:X"))#,note=T,Cnote=T)#,cc("M"))
##
## Model Summary
##
## ───────────────────────────────────────────────────
## (1) Y (2) M
## ───────────────────────────────────────────────────
## (Intercept) 12.589 * 0.940 ***
## (6.046) (0.003)
## Xb -4.048
## (8.692)
## M -9.898
## (6.500)
## W -0.077 0.007 ***
## (0.215) (0.001)
## X -8.992 *
## (4.305)
## Xb:M 3.897
## (9.338)
## Xb:W -0.015
## (0.311)
## W:X 0.180
## (0.122)
## M:X 9.729 *
## (4.611)
## ───────────────────────────────────────────────────
## Marginal R^2 0.085 0.545
## Conditional R^2 0.133 1.000
## AIC 2353.639 -21095.194
## BIC 2411.880 -21077.274
## Num. obs. 652 652
## Num. groups: clus 31 31
## Var: clus (Intercept) 0.244 0.000
## Var: clus X 0.286
## Cov: clus (Intercept) X -0.251
## Var: Residual 2.035 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.025 0.062 3.897 9.338 0.007 0.001 -0.097 0.151 no
## 2 Monte Carlo based Mediation 0.065 0.032 9.729 4.611 0.007 0.001 0.004 0.132 yes
## ────────────────────────────────────────────────────────────────────────────────────
OLS.ab2 <- lm(M ~ W, data = Illu.data, na.action = na.exclude)
#OLS.ab2 <- lm(BA20.localworkerpc~BA20.OSH, data = data, na.action = na.exclude)
#GLM_summary(OLS.ab2)
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) 12.589 * 0.940 ***
## (6.046) (0.004)
## Xb -4.048
## (8.692)
## M -9.898
## (6.500)
## W -0.077 0.006 ***
## (0.215) (0.001)
## X -8.992 *
## (4.305)
## Xb:M 3.897
## (9.338)
## Xb:W -0.015
## (0.311)
## W:X 0.180
## (0.122)
## M:X 9.729 *
## (4.611)
## ────────────────────────────────────────────────
## Marginal R^2 0.085
## Conditional R^2 0.133
## AIC 2353.639
## BIC 2411.880
## Num. obs. 652 652
## Num. groups: clus 31
## Var: clus (Intercept) 0.244
## Var: clus X 0.286
## Cov: clus (Intercept) X -0.251
## Var: Residual 2.035
## R^2 0.038
## Adj. R^2 0.037
## ────────────────────────────────────────────────
## 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.022 0.053 3.897 9.338 0.006 0.001 -0.083 0.130 no
## 2 Monte Carlo based Mediation 0.055 0.028 9.729 4.611 0.006 0.001 0.003 0.115 yes
## ────────────────────────────────────────────────────────────────────────────────────
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 NA
## 2 Monte Carlo based Mediation 0.022 0.053 3.897 9.338 0.006 0.001 -0.083 0.130 no
## 3 Monte Carlo based Mediation 0.055 0.028 9.729 4.611 0.006 0.001 0.003 0.115 yes
## 4 MEDIATION RESULTS BASED ON MLM NA
## 5 Monte Carlo based Mediation 0.025 0.062 3.897 9.338 0.007 0.001 -0.097 0.151 no
## 6 Monte Carlo based Mediation 0.065 0.032 9.729 4.611 0.007 0.001 0.004 0.132 yes
## ───────────────────────────────────────────────────────────────────────────────────────
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.02161 0.05488
## SE 0.05308 0.02811
## ──────────────────────
## Monte Carlo based Estimation (rep=5000, seed=1223 if not preset):
## ────────────────────────────────────────────────────────────────────
## Expression Estimate MCEstimate S.E. LLCI ULCI sig
## ────────────────────────────────────────────────────────────────────
## Monte Carlo p1-p2 -0.033 -0.033 (0.060) -0.150 0.084 no
## ────────────────────────────────────────────────────────────────────
## Parameters and SE:
## ──────────────────────
## p1 p2
## ──────────────────────
## Value 0.02549 0.06463
## SE 0.06212 0.03216
## ──────────────────────
## Monte Carlo based Estimation (rep=5000, seed=1223 if not preset):
## ────────────────────────────────────────────────────────────────────
## Expression Estimate MCEstimate S.E. LLCI ULCI sig
## ────────────────────────────────────────────────────────────────────
## Monte Carlo p1-p2 -0.039 -0.039 (0.070) -0.175 0.098 no
## ────────────────────────────────────────────────────────────────────