1 PREPARATION

1.1 Merge

1.1.1 Analysis 1

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")

1.1.2 Analysis 2

#- 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")

1.2 Loading data

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)

1.3 Model

Theoretical model
Theoretical model
Statistical model
Statistical model

2 ILLUSTRATION-FULL MODEL (PLAN B STUDY 2)

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"))

2.1 MLM-based

2.1.1 ab2 (Mplus approach)

#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)

2.1.2 cw2

#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

2.1.3 Mediation

getRNames(MLM.cw2)
## [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
## ────────────────────────────────────────────────────────────────────────────────────

2.2 OLS-based

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
## ────────────────────────────────────────────────────────────────────────────────────

2.3 Results comparison based on OLS and MLM

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
## ────────────────────────────────────────────────────────────────────
MLM.ce=ESCaculatorMC(params = MLM.Me$ab[1:2], SEs = MLM.Me$SEab[1:2], expression = "p1-p2")
## 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
## ────────────────────────────────────────────────────────────────────

3 MLM INTERACTION ANALYSIS FOR cw1 and cw2

Source: Gardner, R. G., Harris, T. B., Li, N., Kirkman, B. L., & Mathieu, J. E. (2017). Understanding “it depends” in organizational research: A theory-based taxonomy, review, and future research agenda concerning interactive and quadratic relationships. Organizational Research Methods, 20(4), 610-638.
Source: Gardner, R. G., Harris, T. B., Li, N., Kirkman, B. L., & Mathieu, J. E. (2017). Understanding “it depends” in organizational research: A theory-based taxonomy, review, and future research agenda concerning interactive and quadratic relationships. Organizational Research Methods, 20(4), 610-638.

3.1 Plan A

data=dat

###-A. Writing formula
source("https://www.dropbox.com/scl/fi/3hf3qmpsvusakfxzmi8wj/V2.CROSS-LEVEL-INTERACTION-cw.R?rlkey=mg6i0kb7acqfu7vsbmnzls9t3&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="JAPanalysis2a.xlsx", sheetName="Interaction.cw1cw2", append=TRUE, row.names=TRUE)

3.2 Plan B

3.2.1 Study 1

data<- dat[round != 3]

###-A. Writing formula
source("https://www.dropbox.com/scl/fi/3hf3qmpsvusakfxzmi8wj/V2.CROSS-LEVEL-INTERACTION-cw.R?rlkey=mg6i0kb7acqfu7vsbmnzls9t3&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="JAPanalysis2b.xlsx", sheetName="S1Interaction.cw1cw2", append=TRUE, row.names=TRUE)

3.2.2 Study 2

data<- dat[round == 3]

###-A. Writing formula
source("https://www.dropbox.com/scl/fi/3hf3qmpsvusakfxzmi8wj/V2.CROSS-LEVEL-INTERACTION-cw.R?rlkey=mg6i0kb7acqfu7vsbmnzls9t3&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="JAPanalysis2b.xlsx", sheetName="S2Interaction.cw1cw2", append=TRUE, row.names=TRUE)

4 ANALYSIS FOR ab2

Full results: https://www.dropbox.com/scl/fi/4owu44ccfxdliwykmc2na/JAPanalysis2b.xlsx?rlkey=ohjgzicb6k9n9lhwwf8zbfxe6&dl=0

4.1 Plan A

data=dat

##-Running 1-1000
###-A. Writing formula
source("https://www.dropbox.com/scl/fi/ujpghh1p6x2b7bzrajd7n/V2.MAIN-EFFECT-ab2bb2.R?rlkey=d1bgxnc7k8efpr917cxsaq2wt&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="JAPanalysis2a.xlsx", sheetName="Main.ab2bb2", append=TRUE, row.names=TRUE)

4.2 Plan B

Full results: https://www.dropbox.com/scl/fi/4owu44ccfxdliwykmc2na/JAPanalysis2b.xlsx?rlkey=ohjgzicb6k9n9lhwwf8zbfxe6&dl=0

4.2.1 Study 1

data<- dat[round != 3]

##-Running 1-1000
###-A. Writing formula
source("https://www.dropbox.com/scl/fi/ujpghh1p6x2b7bzrajd7n/V2.MAIN-EFFECT-ab2bb2.R?rlkey=d1bgxnc7k8efpr917cxsaq2wt&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="JAPanalysis2b.xlsx", sheetName="S1Main.ab2bb2", append=TRUE, row.names=TRUE)

4.2.2 Study 2

data<- dat[round == 3]

##-Running 1-1000
###-A. Writing formula
source("https://www.dropbox.com/scl/fi/ujpghh1p6x2b7bzrajd7n/V2.MAIN-EFFECT-ab2bb2.R?rlkey=d1bgxnc7k8efpr917cxsaq2wt&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="JAPanalysis2b.xlsx", sheetName="S2Main.ab2bb2", append=TRUE, row.names=TRUE)

5 NO USE

#wY
Wdata=import("C:/Users/Eason Zhang/Dropbox/social impact project/data/employeedata.sav")%>%as.data.table()
#matching_variables.W <- grep(cc("mentalhealth,stress,wage,satoverall"), names(Wdata), value = TRUE)
matching_variables.W <- cc("mentalhealth,stress,wage,satoverall")
contains_variables.W <- names(Wdata) %in% matching_variables.W
matching_variables.W <- names(BA2r.1521)[contains_variables.W]


#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)

#BApractices20192021.sav一共有36个变量。
BA2r.1921=import("C:/Users/Eason Zhang/Dropbox/social impact project/BApractices20192021.sav")%>%as.data.table()
matching_variables.1921 <- BA2r.1921[, grep(paste(targeted_variables, collapse = "|"), 
                               names(BA2r.1921), value = TRUE), with = FALSE]
Nvariables.1921 <- data.table(variable = names(matching_variables.1921))
Nvariables.1921 <- Nvariables.1921[order(variable)]
print_table(Nvariables.1921)