Author

Ana M. Díaz

As a reminder, the two guiding research questions for this second study were:

RQ1: Do parenting dimensions that have been identified in the literature as promotive factors for positive youth outcomes serve the same supportive function for families in informal settlements?

RQ 2: Which individual, relational, and community-level factors are predictors of parental resilience in informal settlements?

1 Preliminary Steps

1.1 Computing and Examining Focal Variables

1.1.1 Item-Level Missingness

Prior to analyses, I established that at least 80% of item-level data needed to be available as the criteria to compute composite variables. [Insert description of how this differed by youth (just 1 obs) and by mother data]. While examining missingness, there was one mother observation with missing values in all items of the ACE questionnaire. Upon further examination, I identified that this participant did not have child-level data (no child in the family was interviewed). I removed this participant, which led to a sample of 100 mothers and 160 youth.

1.1.2 Composite-Level Missingness

Missingness was sparse and occurred predominantly at Level 2 (mother-level variables), ranging from 0.6% to 3.8% across variables.Given the sparse level of missingness at the item level, I followed Huang’s (2022) workflow and steps for conducting a multiple imputation procedure with multilevel data using the jomo, mitml, and naniar packages in R. Prior to imputation, scale scores were treated as missing if more than 20% of constituent items were missing. I first computed scale scores per observation, merged youth and mother datasets by family ID, and applied the multiple imputation procedure to this merged dataset with the study focal variables. The procedure consisted of 20 imputed datasets with 5,000 burn-in iterations to ensure convergence.

I evaluated both SDQ subscales and the Controlling Parenting variable as candidate Level 1 target variables given their intraclass correlation coefficients. The Controlling Parenting variable (ICC = .32) best captured the clustering structure of the data and also yielded the best convergence diagnostics. Convergence was assessed using the Rhat statistic (Huang, 2022), which assesses whether model estimates stabilized over the course of the sampling process, with values above 1.05 indicating non-optimal convergence. All parameters showed adequate convergence (max Rhat ≤ 1.03).

I verified the results of the multiple imputation procedures by creating a data frame with the mean values for each variable across the 20 imputed datasets. I also examined the distribution of these means. As expected, mean values for all Level 1 variables were identical across datasets. There was some variation in Level 2 means across datasets, but they were mostly clustered, suggesting the imputed values were plausible. Variables with higher rates of missingness (biweekly income, 3.8%; maternal education, 1.2%) showed somewhat greater variability, consistent with greater uncertainty in their imputed values.

Code 
#### First creating a function to remove youth data from mom data and viceversa

drop_irrelevant <- function(data, sentinel = c(-77, -88)) {
  data %>% select(where(~ !all(. %in% c(sentinel, NA))))
}


#### Creating a function to assign NA to -77, -88 while keeping the info in 
    # new columns
recode_sentinels <- function(data) {
  data %>%
    mutate(across(everything(), as.vector)) %>%
    mutate(across(
      where(is.numeric),
    list(
      dont_know = ~ if_else(. == -77, 1L, 0L),
      refused = ~ if_else(. == -88, 1L, 0L)
    ),
    .names = "{.fn}__{.col}"
           )) %>%
    mutate(across(
      where(is.numeric), ~ na_if(., -77))) %>%
    mutate(across(
      where(is.numeric), ~ na_if(., -88)))
}


##### Mothers Data #####
mother_data <- subset(data, redcap_repeat_instrument == "" & role == 1) %>%
  drop_irrelevant() %>%
  recode_sentinels()

      ### removing the one observation without child data (record_id = 101, from EDI)
          mother_data <- mother_data %>%
            filter(record_id != 101)

##### Youth Data #####
youth_data <- subset(data, redcap_repeat_instrument >= 1) %>%
  drop_irrelevant() %>%
  recode_sentinels()
          
          #### B. COMPUTING VARIABLES ####
          
rowmeans_threshold <- function(df, threshold = 0.80) {
  apply(df, 1, function(x) {
    prop_present <- mean(!is.na(x))
    if (prop_present >= threshold) mean(x, na.rm = TRUE) else NA_real_
  })
}

##### 1. Mother Mental Health ####
mother_data <- mother_data %>%
  mutate(mother_dsmtotal = rowmeans_threshold(select(., starts_with("dsm_"), -dsm_sud3)))

psych::describe(mother_data$mother_dsmtotal)

##### 2. Mother Social Connections ####
mother_data <- mother_data %>%
  mutate(mother_sc = rowmeans_threshold(select(., starts_with("papf_sc")))) 

psych::describe(mother_data$mother_sc)

##### 3. Community Mutual Aid ####
mother_data <- mother_data %>%
  mutate(mutual_aid = rowmeans_threshold(
    select(., vecinos_apoyo1:vecinos_apoyo7)  
  ))

psych::describe(mother_data$mutual_aid)

##### 4. Comm Participation and Org ####
mother_data <- mother_data %>%
  mutate(
    vecinos_org5_r = 4 - vecinos_org5,
    vecinos_org6_r = 4 - vecinos_org6
  )

mother_data <- mother_data %>%
  mutate(commparticipation = rowmeans_threshold(
    select(., vecinos_org1, vecinos_org2, vecinos_org5_r,
             vecinos_org3, vecinos_org4, vecinos_org6_r, vecinos_org8)
  ))

psych::describe(mother_data$commparticipation)

##### 5. Youth-Reported Caregiving Practices ####

###### 5.1. comunicativo #####
youth_data <- youth_data %>%
  mutate(crpbi_comunicativo = rowmeans_threshold(
    select(., c(crpbi_5, crpbi_6, crpbi_7, crpbi_8, crpbi_9, crpbi_10, crpbi_11)) 
  ))

###### 5.2. hostil/rechazo ####
youth_data <- youth_data %>%
  mutate(crpbi_hostil_rechazo = rowmeans_threshold(
    select(., c(crpbi_21, crpbi_22, crpbi_23, crpbi_24, crpbi_25))  
  ))

###### 5.3. controlador ####
youth_data <- youth_data %>%
  mutate(crpbi_controlador = rowmeans_threshold(
    select(., c(crpbi_16, crpbi_17, crpbi_18, crpbi_19, crpbi_20))
  ))

##### 6. Youth Internalizing and Externalizing Symptoms ####
youth_data <- youth_data %>%
  mutate(
    child_sdq7_r  = 2 - child_sdq7,
    child_sdq21_r = 2 - child_sdq21,
    child_sdq25_r = 2 - child_sdq25,
    child_sdq11_r = 2 - child_sdq11,
    child_sdq14_r = 2 - child_sdq14
  )
          youth_data <- youth_data %>%
            mutate(
              # Emotional problems (3,8,13,16,24)
              sdq_emotional = rowSums(select(., c(child_sdq3, child_sdq8, child_sdq13, child_sdq16, child_sdq24)), na.rm = TRUE),
              
              # Conduct problems (5,7[r],12,18,22)
              sdq_conduct = rowSums(select(., c(child_sdq5, child_sdq7_r, child_sdq12, child_sdq18, child_sdq22)), na.rm = TRUE),
              
              # Hyperactivity (2,10,15,21[r],25[r])
              sdq_hyperactivity = rowSums(select(., c(child_sdq2, child_sdq10, child_sdq15, child_sdq21_r, child_sdq25_r)), na.rm = TRUE),
              
              # Peer problems (6,11[r],14[r],19,23)
              sdq_peer = rowSums(select(., c(child_sdq6, child_sdq11_r, child_sdq14_r, child_sdq19, child_sdq23)), na.rm = TRUE),
              
              # Prosocial (1,4,9,17,20)
              sdq_prosocial = rowSums(select(., c(child_sdq1, child_sdq4, child_sdq9, child_sdq17, child_sdq20)), na.rm = TRUE),
              
              # Additional composite scores
              sdq_externalising = sdq_conduct + sdq_hyperactivity,
              sdq_internalising = sdq_emotional + sdq_peer,
              
              # Total difficulties (excluding prosocial)
              sdq_total = sdq_emotional + sdq_conduct + sdq_hyperactivity + sdq_peer
            )
          
          psych::describe(youth_data$sdq_total)
          
          ##### 7. Mom ACES  ####
          ## dichotomizing
          mother_data <- mother_data %>%
            mutate(across(c(aces_3,
                            aces_7,
                            aces_9,
                            aces_19,
                            aces_20,
                            aces_21,
                            aces_22,
                            aces_23,
                            aces_24,
                            aces_25
            ), 
            ~ as.numeric(. > 1), 
            .names = "{.col}_dummy"))
          
          mother_data <- mother_data %>%
            mutate(aces_1r = 3 - aces_1)
          
          
          
          ### computing index of total aces
          mother_data <- mother_data %>%
            mutate(mother_aces = rowSums(select(., c(
              aces_1r,
              aces_3_dummy,
              aces_7_dummy,
              aces_9_dummy,
              aces_19_dummy,
              aces_20_dummy,
              aces_21_dummy,
              aces_22_dummy,
              aces_23_dummy,
              aces_24_dummy,
              aces_25_dummy,
              aces_10,
              aces_12,
              aces_13,
              aces_14,
              aces_16,
              aces_17,
              aces_18,
              aces_26,
              aces_27,
              aces_28
            )), na.rm = TRUE))
          
          
          
          mother_data <- mother_data %>%
            mutate(
              # Create domain scores
              ace_sexual_abuse = ifelse(
                aces_26 == 1 | aces_27 == 1 | aces_28 == 1, 1, 0),
              
              ace_physical_abuse = ifelse(
                aces_21_dummy == 1 | aces_23_dummy == 1 | aces_24_dummy == 1, 1, 0),
              
              ace_emotional_abuse = ifelse(
                aces_19_dummy == 1 | aces_20_dummy == 1 | 
                  aces_22_dummy == 1 | aces_25_dummy == 1, 1, 0)
            ) %>%
            rowwise() %>%
            mutate(
              mother_aces = sum(c(ace_sexual_abuse, ace_physical_abuse, ace_emotional_abuse,
                                  aces_10, aces_12, aces_13, aces_14, aces_16, aces_17, aces_18,
                                  aces_1r, aces_3_dummy, aces_7_dummy, aces_9_dummy), 
                                na.rm = TRUE)
            ) %>%
            ungroup()
          
          psych::describe(mother_data$mother_aces)
          
          ##### 8. IS INDEX  ####
          mother_data <- mother_data %>%
            mutate(
              ###### Infrastructure
              water_insecurity = ifelse(agua_acceso_frecuencia %in% c(4, 5, 6), 1, 0),
              waste_inadequate = ifelse(basura_manejo %in% c(1, 2), 1, 0),
              electricity_cutoff = ifelse(electricidad_cortes1 == 1, 1, 0),
              hygiene_dangers = ifelse(
                higienicos_eventos1 >= 1 | higienicos_eventos2 >= 1 | 
                  higienicos_eventos3 >= 1 | higienicos_eventos4 >= 1 | 
                  higienicos_eventos5 >= 1, 1, 0),
              
              ###### Environmental vulnerabilities
              flooding = ifelse(inundaciones > 0, 1, 0),
              
              ###### Tenure insecurity
              no_land_title = ifelse(terrenos2 == 0, 1, 0),
              eviction_exp = ifelse(terrenos4 == 1, 1, 0),
              
              ###### Community safety
              community_violence = ifelse(
                seg_eventos4 >= 1 | seg_eventos5 >= 1 | seg_eventos11 >= 1 |
                  seg_eventos9 >= 1 | seg_eventos12 >= 1 | seg_eventos13 >= 1 | 
                  seg_eventos14 >= 1, 1, 0)
            ) %>%
            ###### Create the overall index using available data
          rowwise() %>%
            mutate(
              informal_settlement_index = sum(c(water_insecurity, waste_inadequate, 
                                                electricity_cutoff, flooding, hygiene_dangers, no_land_title, 
                                                eviction_exp, community_violence), na.rm = TRUE)
            ) %>%
            ungroup()
          
          psych::describe(mother_data$informal_settlement_index)
          
##### creating data sets to merge for imputation #####
youth_data_focal <- youth_data %>%
  select(record_id, sdq_externalising, sdq_internalising, sdq_total,
         crpbi_comunicativo, crpbi_controlador, crpbi_hostil_rechazo,
         child_demo_age,
         child_demo_gender)

youth_data_focal <- youth_data_focal %>%
  mutate(child_demo_gender = factor(child_demo_gender, ordered = FALSE))

mother_data_focal <- mother_data %>%
  select(record_id, mother_aces, informal_settlement_index,
         mother_dsmtotal, mother_sc, mutual_aid, commparticipation,
         adult_education, quincena
  )

# Join into one long dataset where mother variable repeats for each child
merged_focal <- youth_data_focal %>%
  left_join(mother_data_focal, by = "record_id")

# Check structure looks right
merged_focal %>% 
  count(record_id) %>% 
  summary() 
Percent missing per variable:
                record_id         sdq_externalising         sdq_internalising 
                  0.00000                   0.00000                   0.00000 
                sdq_total        crpbi_comunicativo         crpbi_controlador 
                  0.00000                   0.01250                   0.00000 
     crpbi_hostil_rechazo            child_demo_age         child_demo_gender 
                  0.00625                   0.00000                   0.00000 
              mother_aces informal_settlement_index           mother_dsmtotal 
                  0.00625                   0.00625                   0.00625 
                mother_sc                mutual_aid         commparticipation 
                  0.00625                   0.00625                   0.00625 
          adult_education                  quincena 
                  0.01250                   0.03750 

Percent complete cases: 0.94375 
(Minimum) number to impute: 6

Figure 1: Matrix of Missing Values

Table X
Descriptive Statistics for Observed Study Variables
Variable n M SD Min Max
Child Demographics
Child Age 160 11.88 2.95 7.00 19.00
Child Outcomes
SDQ Externalising 160 7.05 3.52 0.00 16.00
SDQ Internalising 160 6.72 3.43 0.00 16.00
SDQ Total 160 13.77 5.46 3.00 30.00
Parenting
Communicative Parenting 158 2.63 0.34 1.29 3.00
Controlling Parenting 160 2.14 0.49 1.00 3.00
Hostile/Rejecting Parenting 159 1.62 0.43 1.00 3.00
Mother Characteristics
Maternal Education 98 9.83 3.51 1.00 15.00
Informal Settlement Index 99 4.76 1.21 2.00 8.00
Maternal ACEs 99 4.20 3.36 0.00 14.00
DSM Cross-Cutting Symptom 99 0.83 0.69 0.00 2.91
Mother Social Connections 99 2.90 1.03 0.00 4.00
Biweekly Income 96 316.84 152.50 0.00 800.00
Community
Community Participation 99 2.62 0.62 0.57 4.00
Perceived Mutual Aid 99 2.78 0.64 1.14 4.00
Note. Statistics for child-level variables (Child Outcomes, Parenting, Child Demographics) are based on the full sample of youth observations. Statistics for mother-level variables (Mother Characteristics, Community) are based on one observation per family.

1.2 Multiple Imputation

To assess the results (convergence and efficiency) of the multiple imputation procedure, I relied on the Rhat statistic, again following Huang (2022), which assesses whether the model’s estimates stabilized over the course of the sampling process (with values above 1.05 indicating non-optimal convergence). All estimates showed adequate convergence (Rhat < 1.05), with the exception of the Level 1 random effects covariance (Psi; max Rhat = 1.21)

I verified the results of the multiple imputation procedures by examining the distribution of variable means across the 20 imputed datasets. As expected, mean values for all Level 1 variables were identical across datasets, reflecting the absence of missing data at that level. There was some variation in Level 2 means across datasets, but values were tightly clustered, suggesting stable and plausible imputed values. Variables with higher rates of missingness (biweekly income, 3.8%; maternal education, 1.2%) showed somewhat greater variability across datasets, consistent with greater uncertainty in their imputed values.

Code 
fml <- list(
  crpbi_controlador + crpbi_comunicativo + crpbi_hostil_rechazo ~ 
    1 + sdq_externalising + sdq_internalising + 
    child_demo_age + child_demo_gender + (1|record_id),
  
  mother_aces + informal_settlement_index + mother_dsmtotal + 
    mother_sc + mutual_aid + commparticipation + 
    adult_education + quincena ~ 1
)

imp <- jomoImpute(
  data = merged_focal,
  formula = fml,
  m = 20,
  seed = 123,
  n.burn = 5000
)

summary(imp)

implist <- mitmlComplete(imp, print = "all")

implist <- within(implist, {
  child_demo_gender <- relevel(as.factor(child_demo_gender), ref = "0")
  sdq_ext_gpm <- clusterMeans(sdq_externalising, record_id)
  sdq_int_gpm <- clusterMeans(sdq_internalising, record_id)
})

# Getting means from each dataset and summaring across them
means_list <- lapply(implist, function(d) {
  colMeans(d[, c("mother_aces", "informal_settlement_index", 
                 "mother_dsmtotal", "mother_sc", "mutual_aid", 
                 "commparticipation", "adult_education", "quincena")], 
           na.rm = TRUE)
})

# Convert to dataframe — rows = datasets, columns = variables
means_df <- do.call(rbind, means_list)
rownames(means_df) <- paste0("imp_", 1:20)

# Summary across all 20 imputations
round(means_df, 3)

# And the overall average across all 20
round(colMeans(means_df), 3)

### Checking youth variables
youth_vars <- c("sdq_externalising", "sdq_internalising", 
                "crpbi_comunicativo", "crpbi_controlador", "crpbi_hostil_rechazo",
                "child_demo_age")

youth_means <- do.call(rbind, lapply(implist, function(d) {
  colMeans(d[, youth_vars], na.rm = TRUE)
}))

rownames(youth_means) <- paste0("imp_", 1:20)
round(youth_means, 3)
Code 
# Overall mean across the 20 imputed datasets
imputed_means <- colMeans(means_df)
imputed_means_df <- data.frame(
  variable = names(imputed_means),
  imp_mean = as.numeric(imputed_means)
)

# Pivot to long
means_long <- as.data.frame(means_df) %>%
  mutate(imp = 1:20) %>%
  pivot_longer(-imp, names_to = "variable", values_to = "mean")

# Merge
plot_df <- left_join(means_long, imputed_means_df, by = "variable")

# Label mapping
var_labels <- c(
  mother_aces = "Maternal ACEs",
  informal_settlement_index = "Informal Settlement Index",
  mother_dsmtotal = "DSM Cross-Cutting Symptom",
  mother_sc = "Mother Social Connections",
  mutual_aid = "Perceived Mutual Aid",
  commparticipation = "Community Participation and Organization",
  adult_education = "Maternal Education",
  quincena = "Biweekly Income"
)

# Apply labels to plot_df
plot_df <- plot_df %>%
  mutate(variable = recode(variable, !!!var_labels))

# Update imputed means df too
imputed_means_df <- imputed_means_df %>%
  mutate(variable = recode(variable, !!!var_labels))

# Re-merge
plot_df <- left_join(
  plot_df %>% select(-imp_mean), 
  imputed_means_df, 
  by = "variable"
)

plot_df <- plot_df %>%
  mutate(mean_centered = mean - imp_mean)
Code 
ggplot(plot_df, aes(x = mean_centered)) +
  geom_histogram(bins = 10, fill = "steelblue", color = "white", alpha = 0.8) +
  geom_vline(xintercept = 0, color = "red", linewidth = 0.8, linetype = "dashed") +
  facet_wrap(~ variable, scales = "free", labeller = label_wrap_gen(width = 25)) +
  labs(
    x = "Deviation from average imputed mean",
    y = "Count",
    caption = "Note. Red dashed line = mean across all 20 imputed datasets (centered at 0)"
  ) +
  theme_minimal()

Figure 2: Examining the Stability of Variable Means Across 20 Imputed Datasets

2 RQ 1: Parenting Dimensions and Youth Internalizing and Externalizing Behaviors

Model 0

\[SDQ_{ij} = \beta_{0j} + r_{ij} \tag{Eq. 1}\]

\[\beta_{0j} = \gamma_{00} + u_{0j} \tag{Eq. 2}\]

Model 1

\[SDQ_{ij} = \beta_{0j} + \beta_1(PQuality)_{ij} \tag{Eq. 3}\]

Level 2: Same as Model 0

Model 2

\[SDQ_{ij} = \beta_{0j} + \beta_1(PQuality)_{ij} + \beta_2(ChildAge)_{ij} + \beta_3(ChildFemale)_{ij} + r_{ij} \tag{Eq. 4}\]

Level 2: Same as Model 0

Where:

  • \(SDQ_{ij}\) = Internalizing/Externalizing score for youth \(i\) in family \(j\)
  • \(\beta_1\) = Fixed effect of parent quality dimension
  • \(\beta_2\) = Fixed effect of youth age
  • \(\beta_3\) = Fixed effect of youth sex (female)
  • \(ChildAge_{ij}\) = Age of youth \(i\) in family \(j\)
  • \(ChildFemale_{ij}\) = 1 if youth \(i\) in family \(j\) is female
  • \(r_{ij}\) = Level 1 residual (individual youth deviation from family mean)
  • \(\gamma_{00}\) = Fixed intercept (grand mean of internalizing/externalizing across all youth)
  • \(u_{0j}\) = Level 2 residual (random family intercept; family \(j\)’s deviation from the grand mean)

Figure 3: Scatterplot of Relationship between Parenting Dimensions and Internalizing

Figure 4: Scatterplot of Relationship between Parenting Dimensions and Externalizing

I then used multilevel regression to assess bivariate associations between youth psychosocial health and parenting dimensions. Specifically, I modeled internalizing and externalizing symptoms on three parenting dimensions (warmth/communicative, hostility, control). After examining these bivariate associations, I included youth age and youth sex in these models.

Code 
###### RQ1: ICCs ######
# Null models - no predictors
null_model_int <- with(implist_std, lmer(sdq_internalising ~ 1 + (1|record_id)))
null_model_ext <- with(implist_std, lmer(sdq_externalising ~ 1 + (1|record_id)))

testEstimates(null_model_ext, extra.pars = TRUE)
testEstimates(null_model_int, extra.pars = TRUE)

###### RQ1.1.EXTERNALIZING #####
    model_int_comunicativo1 <- with(implist_std, lmer(sdq_internalising ~ crpbi_comunicativo + (1 | record_id), REML = FALSE))
    
    model_int_comunicativo2 <- with(implist_std, lmer(sdq_internalising ~ crpbi_comunicativo + 
                                 child_demo_age + child_demo_gender +
                                 (1 | record_id), REML = FALSE))
    
        #### Int. Comunicativo: Pooling estimates ####
                testEstimates(model_int_comunicativo1, extra.pars = TRUE)
                testEstimates(model_int_comunicativo2, extra.pars = TRUE)

    
    model_int_controlador1 <- with(implist_std, lmer(sdq_internalising ~ crpbi_controlador + (1 | record_id), 
                              REML = FALSE))
    model_int_controlador2  <- with(implist_std, lmer(sdq_internalising ~ crpbi_controlador + 
                                     child_demo_age + child_demo_gender +
                                     (1 | record_id), REML = FALSE))
    #### Int. Controlador: Pooling estimates ####    
          testEstimates(model_int_controlador1, extra.pars = TRUE)
          testEstimates(model_int_controlador2, extra.pars = TRUE)

    model_int_rechazo1 <- with(implist_std, lmer(sdq_internalising ~ crpbi_hostil_rechazo + (1 | record_id), 
                          REML = FALSE))
    model_int_rechazo2 <- with(implist_std, lmer(sdq_internalising ~ crpbi_hostil_rechazo + 
                                     child_demo_age + child_demo_gender +
                                     (1 | record_id), REML = FALSE))
    
    #### Int. Rechazo: Pooling estimates ####    
            testEstimates(model_int_rechazo1, extra.pars = TRUE)
            testEstimates(model_int_rechazo2, extra.pars = TRUE)
   
    
    
###### RQ1.2.EXTERNALIZING ######
    model_ext_comunicativo1 <- with(implist_std, lmer(sdq_externalising ~ crpbi_comunicativo + (1 | record_id), 
                               REML = FALSE))
    
    model_ext_comunicativo2 <- with(implist_std, lmer(sdq_externalising ~ crpbi_comunicativo + 
                                     child_demo_age + child_demo_gender +
                                     (1 | record_id), REML = FALSE))
    
            #### Ext. Comunicativo: Pooling estimates ####
            testEstimates(model_ext_comunicativo1, extra.pars = TRUE)
            testEstimates(model_ext_comunicativo2, extra.pars = TRUE)
    
    model_ext_controlador1 <- with(implist_std, lmer(sdq_externalising ~ crpbi_controlador + (1 | record_id), 
                           REML = FALSE))
    
    model_ext_controlador2  <- with(implist_std, lmer(sdq_externalising ~ crpbi_controlador + 
                                     child_demo_age + child_demo_gender +
                                     (1 | record_id), REML = FALSE))
    
                  #### Extt. Controlador: Pooling estimates ####    
                  testEstimates(model_ext_controlador1, extra.pars = TRUE)
                  testEstimates(model_ext_controlador2, extra.pars = TRUE)

    model_ext_rechazo1 <- with(implist_std, lmer(sdq_externalising ~ crpbi_hostil_rechazo + (1 | record_id), 
                          REML = FALSE))
    
    model_ext_rechazo2 <- with(implist_std, lmer(sdq_externalising ~ crpbi_hostil_rechazo + 
                                 child_demo_age + child_demo_gender +
                                 (1 | record_id), REML = FALSE))
    
    #### Ext. Rechazo: Pooling estimates ####    
          testEstimates(model_ext_rechazo1, extra.pars = TRUE)
          testEstimates(model_ext_rechazo2, extra.pars = TRUE)

[I need to insert comparison of the standardized regression coefficients in the bivariate models to effect sizes in the Pinquart metanalyses.]

?(caption)

Table 9
Pooled Multilevel Model Estimates: Communicative Parenting
Model 1
Model 2
Model 1
Model 2
Predictor
Internalising Behaviour
Externalising Behaviour
≈ [95% CI] SE p β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p
Fixed Effects
Intercept 0.007 [-0.161, 0.175] 0.086 0.934 -0.176 [-0.387, 0.035] 0.108 0.102 -0.005 [-0.164, 0.154] 0.081 0.950 0.121 [-0.090, 0.333] 0.108 0.260
Communicative Parenting -0.079 [-0.232, 0.075] 0.078 0.314 -0.092 [-0.238, 0.055] 0.075 0.219 -0.169 [-0.322, -0.015] 0.078 0.031 -0.171 [-0.326, -0.017] 0.079 0.030
Child Age — — — -0.229 [-0.375, -0.083] 0.074 0.002 — — — 0.051 [-0.102, 0.203] 0.078 0.514
Child Sex — — — 0.382 [0.098, 0.665] 0.145 0.008 — — — -0.268 [-0.568, 0.032] 0.153 0.079
Variance Components
τ00 (between-cluster) 0.229 — — 0.253 — — 0.104 — — 0.103 — —
σ² (within-cluster) 0.762 — — 0.651 — — 0.859 — — 0.838 — —
ICC 0.231 — — 0.280 — — 0.108 — — 0.110 — —
Note. β = standardized coefficient. 95% CI computed as β ± 1.96 × SE. Estimates pooled across 20 imputed datasets. τ00 = random intercept variance. σ² = residual variance. ICC = intraclass correlation coefficient. Bold values indicate statistical significance (p < .05).
Table 10
Pooled Multilevel Model Estimates: Controlling Parenting
Model 1
Model 2
Model 1
Model 2
Predictor
Internalising Behaviour
Externalising Behaviour
≈ [95% CI] SE p β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p
Fixed Effects
Intercept 0.007 [-0.158, 0.171] 0.084 0.937 -0.279 [-0.481, -0.078] 0.103 0.007 -0.008 [-0.168, 0.152] 0.082 0.922 0.057 [-0.158, 0.272] 0.110 0.606
Controlling Parenting 0.325 [0.178, 0.472] 0.075 < .001 0.384 [0.236, 0.532] 0.076 < .001 0.182 [0.029, 0.336] 0.078 0.020 0.191 [0.028, 0.353] 0.083 0.021
Child Age — — — -0.129 [-0.267, 0.010] 0.071 0.069 — — — 0.121 [-0.035, 0.276] 0.079 0.129
Child Sex — — — 0.601 [0.328, 0.874] 0.139 < .001 — — — -0.133 [-0.445, 0.179] 0.159 0.403
Variance Components
τ00 (between-cluster) 0.274 — — 0.267 — — 0.122 — — 0.106 — —
σ² (within-cluster) 0.636 — — 0.532 — — 0.839 — — 0.834 — —
ICC 0.301 — — 0.334 — — 0.127 — — 0.112 — —
Note. β = standardized coefficient. 95% CI computed as β ± 1.96 × SE. Estimates pooled across 20 imputed datasets. τ00 = random intercept variance. σ² = residual variance. ICC = intraclass correlation coefficient. Bold values indicate statistical significance (p < .05).
Table 11
Pooled Multilevel Model Estimates: Hostile/Rejecting Parenting
Model 1
Model 2
Model 1
Model 2
Predictor
Internalising Behaviour
Externalising Behaviour
≈ [95% CI] SE p β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p
Fixed Effects
Intercept 0.004 [-0.158, 0.167] 0.083 0.957 -0.214 [-0.418, -0.011] 0.104 0.039 -0.004 [-0.149, 0.141] 0.074 0.956 0.067 [-0.127, 0.262] 0.099 0.499
Hostile/Rejecting Parenting 0.248 [0.100, 0.397] 0.076 0.001 0.256 [0.115, 0.398] 0.072 < .001 0.393 [0.250, 0.535] 0.073 < .001 0.395 [0.253, 0.537] 0.073 < .001
Child Age — — — -0.192 [-0.332, -0.051] 0.072 0.007 — — — 0.114 [-0.028, 0.255] 0.072 0.114
Child Sex — — — 0.459 [0.183, 0.734] 0.140 0.001 — — — -0.149 [-0.431, 0.134] 0.144 0.303
Variance Components
τ00 (between-cluster) 0.206 — — 0.232 — — 0.047 — — 0.033 — —
σ² (within-cluster) 0.726 — — 0.612 — — 0.789 — — 0.783 — —
ICC 0.221 — — 0.275 — — 0.056 — — 0.040 — —
Note. β = standardized coefficient. 95% CI computed as β ± 1.96 × SE. Estimates pooled across 20 imputed datasets. τ00 = random intercept variance. σ² = residual variance. ICC = intraclass correlation coefficient. Bold values indicate statistical significance (p < .05).

3 RQ2: Factors that Promote Parental Resilience

3.1 Regression Assumptions

[finalize regression assumption assessments]

Prior to running the conditional models, I generated a correlation matrix between all predictors for the final model. The goal was to assess potential collinearity issues as well as to examine the qualities of the relationship between focal study variables.

Code 
# Select parenting composites and main predictors
cor_vars <- merged_focal %>%
  select(crpbi_comunicativo, 
         crpbi_controlador, 
         crpbi_hostil_rechazo,
         mother_aces, 
         informal_settlement_index, 
         mother_dsmtotal,
         mother_sc, 
         commparticipation, 
         mutual_aid,
         adult_education,
         child_demo_age,
         sdq_total)

# Basic correlation matrix
cor_matrix <- cor(cor_vars, use = "pairwise.complete.obs")
round(cor_matrix, 2)

# With p-values
library(Hmisc)
rcorr(as.matrix(cor_vars))

# For a nicer visualization
library(corrplot)
corrplot(cor_matrix, method = "color", type = "upper",
         addCoef.col = "black", number.cex = 0.7,
         tl.cex = 0.8, tl.col = "black")

I also examined the eigenvalues of the level 2 predictor-only correlation matrix. I previously determined that if the sum of their reciprocal was larger than 35 (five times seven predictors) this would indicate issues of collinearity (Zieffler, 2023).

[I still need to extract and insert eigenvaleus matrix analyses]

3.2 Multilevel Model Specification

Following these assessments of model assumptions, four sequential models tested unique variance explained by each ecological level. The goal with this approach was to examine whether more distal factors (childhood adversity, structural conditions, community resources) explained variance in parenting quality beyond more proximal factors (current mental health, relational resources), while accounting for youth characteristics. Given that these are multilevel models, I needed to account for both within family (youth level) and between family variance (mother-level). All these models were fitted using the 20 imputed datasets.

  • Model 1 included youth-level covariates (age, sex, total mental health problems) and family demographics (maternal education).
  • Model 2 added mother-level factors (mental health, ACEs, informal settlement adversity).
  • Model 3 added relational-level factors (social connections).
  • Model 4 added community-level factors (mutual aid, community participation).

All Level 2 predictors were grand mean centered and the models were estimated employing Maximum Likelihood Estimation (MLE) using the lme4 package in R (Bates et al., 2015) to allow for model comparisons (likelihood ratio tests). I then estimated the final models with reduced maximum likelihood (REML) due to its appropriateness for estimating random components in multilevel models (Finch et al., 2019). I used these final estimates to compare with the standard OLS estimates in the final step of analyses.

Code 
###### RQ.1. Comunicativo #####
        com_m0 <- with(implist_std, lmer(crpbi_comunicativo ~ 1 + (1 | record_id),  REML = FALSE))
        
        com_m1 <- with(implist_std, lmer(crpbi_comunicativo ~ sdq_total + child_demo_age + child_demo_gender +
                         adult_education +
                         (1 | record_id), REML = FALSE))
        
        com_m2 <- with(implist_std, lmer(crpbi_comunicativo ~ sdq_total + child_demo_age + child_demo_gender +
                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                         (1 | record_id), REML = FALSE))
        
        com_m3 <- with(implist_std, lmer(crpbi_comunicativo ~ sdq_total + child_demo_age + child_demo_gender +
                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                         mother_sc +
                         (1 | record_id), REML = FALSE))
        
        com_m4_ml <- with(implist_std, lmer(crpbi_comunicativo ~ sdq_total + child_demo_age + child_demo_gender +
                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                         mother_sc + commparticipation + mutual_aid +
                         (1 | record_id), REML = FALSE))

                #### Comunicativo - Pooling estimates
                          testEstimates(com_m1, extra.pars = TRUE)
                          testEstimates(com_m2, extra.pars = TRUE)
                          testEstimates(com_m3, extra.pars = TRUE)
                          testEstimates(com_m4_ml, extra.pars = TRUE)

###### RQ.2. Controlador #####
          ctrl_m0 <- with(implist_std, lmer(crpbi_controlador ~ 1 + (1 | record_id),  REML = FALSE))
          
          ctrl_m1 <- with(implist_std, lmer(crpbi_controlador ~ sdq_total + child_demo_age + child_demo_gender +
                                         adult_education +
                                         (1 | record_id), REML = FALSE))
          
          ctrl_m2 <- with(implist_std, lmer(crpbi_controlador ~ sdq_total + child_demo_age + child_demo_gender +
                                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                         (1 | record_id), REML = FALSE))
          
          ctrl_m3 <- with(implist_std, lmer(crpbi_controlador ~ sdq_total + child_demo_age + child_demo_gender +
                                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                         mother_sc +
                                         (1 | record_id), REML = FALSE))
          
          ctrl_m4_ml <- with(implist_std, lmer(crpbi_controlador ~ sdq_total + child_demo_age + child_demo_gender +
                                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                         mother_sc + commparticipation + mutual_aid +
                                         (1 | record_id), REML = FALSE))
          
          #### Controlador - Pooling estimates
          testEstimates(ctrl_m1, extra.pars = TRUE)
          testEstimates(ctrl_m2, extra.pars = TRUE)
          testEstimates(ctrl_m3, extra.pars = TRUE)
          testEstimates(ctrl_m4_ml, extra.pars = TRUE)


###### RQ.3. Hostil/Rechazo #####
          hos_m0 <- with(implist_std, lmer(crpbi_hostil_rechazo ~ 1 + (1 | record_id),  REML = FALSE))
          
          hos_m1 <- with(implist_std, lmer(crpbi_hostil_rechazo ~ sdq_total + child_demo_age + child_demo_gender +
                                          adult_education +
                                          (1 | record_id), REML = FALSE))
          
          hos_m2 <- with(implist_std, lmer(crpbi_hostil_rechazo ~ sdq_total + child_demo_age + child_demo_gender +
                                          adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                          (1 | record_id), REML = FALSE))
          
          hos_m3 <- with(implist_std, lmer(crpbi_hostil_rechazo ~ sdq_total + child_demo_age + child_demo_gender +
                                          adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                          mother_sc +
                                          (1 | record_id), REML = FALSE))
          
          hos_m4_ml <- with(implist_std, lmer(crpbi_hostil_rechazo ~ sdq_total + child_demo_age + child_demo_gender +
                                          adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                          mother_sc + commparticipation + mutual_aid +
                                          (1 | record_id), REML = FALSE))
          
          #### Hostil/Rechazo - Pooling estimates
          testEstimates(hos_m1, extra.pars = TRUE)
          testEstimates(hos_m2, extra.pars = TRUE)
          testEstimates(hos_m3, extra.pars = TRUE)
          testEstimates(hos_m4_ml, extra.pars = TRUE)
Code 
    com_m4_reml <- with(implist_std, lmer(crpbi_comunicativo ~ sdq_total + child_demo_age + child_demo_gender +
                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                         mother_sc + commparticipation + mutual_aid +
                         (1 | record_id), REML = TRUE))

     ctrl_m4_reml <- with(implist_std, lmer(crpbi_controlador ~ sdq_total + child_demo_age + child_demo_gender +
                                         adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                         mother_sc + commparticipation + mutual_aid +
                                         (1 | record_id), REML = TRUE))
     
     hos_m4_reml <- with(implist_std, lmer(crpbi_hostil_rechazo ~ sdq_total + child_demo_age + child_demo_gender +
                                          adult_education + mother_aces + informal_settlement_index + mother_dsmtotal +
                                          mother_sc + commparticipation + mutual_aid +
                                          (1 | record_id), REML = TRUE))
Table 1
Pooled Multilevel Model Estimates for Communicative Parenting
Predictor
Model 1
Model 2
Model 3
Model 4
β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p
Fixed Effects
Intercept 0.063 [-0.145, 0.270] 0.106 0.554 0.051 [-0.157, 0.259] 0.106 0.631 0.050 [-0.158, 0.258] 0.106 0.636 0.048 [-0.160, 0.256] 0.106 0.649
SDQ Total -0.174 [-0.330, -0.019] 0.079 0.028 -0.184 [-0.342, -0.025] 0.081 0.023 -0.183 [-0.342, -0.025] 0.081 0.023 -0.188 [-0.347, -0.029] 0.081 0.021
Child Age -0.156 [-0.310, -0.002] 0.079 0.048 -0.146 [-0.300, 0.008] 0.079 0.063 -0.146 [-0.300, 0.008] 0.079 0.063 -0.149 [-0.303, 0.006] 0.079 0.060
Child Sex -0.133 [-0.437, 0.172] 0.155 0.393 -0.108 [-0.416, 0.200] 0.157 0.491 -0.107 [-0.415, 0.202] 0.157 0.498 -0.103 [-0.412, 0.207] 0.158 0.514
Maternal Education 0.037 [-0.116, 0.190] 0.078 0.632 0.057 [-0.101, 0.214] 0.080 0.480 0.063 [-0.100, 0.225] 0.083 0.450 0.050 [-0.116, 0.216] 0.085 0.552
Maternal ACEs — — — 0.070 [-0.097, 0.237] 0.085 0.410 0.074 [-0.095, 0.244] 0.086 0.390 0.078 [-0.092, 0.247] 0.087 0.371
Informal Settlement Index — — — -0.055 [-0.216, 0.106] 0.082 0.502 -0.054 [-0.215, 0.107] 0.082 0.511 -0.060 [-0.222, 0.101] 0.083 0.465
DSM Cross-Cutting Symptom — — — -0.086 [-0.254, 0.082] 0.086 0.316 -0.082 [-0.252, 0.088] 0.087 0.346 -0.080 [-0.250, 0.090] 0.087 0.359
Mother Social Connections — — — — — — 0.024 [-0.141, 0.190] 0.084 0.773 0.031 [-0.138, 0.200] 0.086 0.718
Community Participation — — — — — — — — — 0.033 [-0.122, 0.188] 0.079 0.676
Perceived Mutual Aid — — — — — — — — — -0.050 [-0.216, 0.115] 0.085 0.552
Variance Components
τ00 (between-cluster) 0.005 — — 0.003 — — 0.003 — — 0.001 — —
σ² (within-cluster) 0.932 — — 0.924 — — 0.924 — — 0.923 — —
ICC 0.005 — — 0.003 — — 0.003 — — 0.001 — —
Note. β = Standardized coefficient. 95% CI computed as β ± 1.96 × SE. Estimates pooled across 20 imputed datasets. τ00 = random intercept variance. σ² = residual variance. ICC = intraclass correlation coefficient. Bold values indicate statistical significance (p < .05).
Table 2
Pooled Multilevel Model Estimates for Controlling Parenting
Predictor
Model 1
Model 2
Model 3
Model 4
β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p
Fixed Effects
Intercept 0.270 [0.074, 0.466] 0.100 0.007 0.273 [0.078, 0.467] 0.099 0.006 0.272 [0.077, 0.466] 0.099 0.006 0.275 [0.080, 0.469] 0.099 0.006
SDQ Total 0.324 [0.189, 0.459] 0.069 < .001 0.324 [0.188, 0.460] 0.069 < .001 0.324 [0.188, 0.460] 0.069 < .001 0.323 [0.187, 0.459] 0.070 < .001
Child Age -0.213 [-0.345, -0.081] 0.067 0.002 -0.211 [-0.342, -0.080] 0.067 0.002 -0.211 [-0.342, -0.080] 0.067 0.002 -0.214 [-0.346, -0.083] 0.067 0.001
Child Sex -0.564 [-0.817, -0.310] 0.129 < .001 -0.565 [-0.818, -0.311] 0.129 < .001 -0.563 [-0.816, -0.310] 0.129 < .001 -0.567 [-0.821, -0.312] 0.130 < .001
Maternal Education -0.080 [-0.239, 0.078] 0.081 0.321 -0.100 [-0.259, 0.059] 0.081 0.216 -0.094 [-0.258, 0.071] 0.084 0.265 -0.104 [-0.272, 0.064] 0.086 0.226
Maternal ACEs — — — 0.104 [-0.067, 0.275] 0.087 0.232 0.109 [-0.065, 0.284] 0.089 0.218 0.110 [-0.064, 0.283] 0.089 0.215
Informal Settlement Index — — — 0.044 [-0.119, 0.206] 0.083 0.598 0.044 [-0.118, 0.206] 0.083 0.595 0.039 [-0.124, 0.201] 0.083 0.641
DSM Cross-Cutting Symptom — — — 0.069 [-0.108, 0.245] 0.090 0.445 0.073 [-0.105, 0.251] 0.091 0.422 0.075 [-0.103, 0.252] 0.091 0.408
Mother Social Connections — — — — — — 0.027 [-0.142, 0.195] 0.086 0.756 0.027 [-0.143, 0.198] 0.087 0.752
Community Participation — — — — — — — — — 0.052 [-0.105, 0.210] 0.080 0.517
Perceived Mutual Aid — — — — — — — — — -0.032 [-0.198, 0.134] 0.085 0.705
Variance Components
τ00 (between-cluster) 0.298 — — 0.272 — — 0.273 — — 0.266 — —
σ² (within-cluster) 0.482 — — 0.482 — — 0.481 — — 0.483 — —
ICC 0.382 — — 0.361 — — 0.362 — — 0.355 — —
Note. β = Standardized coefficient. 95% CI computed as β ± 1.96 × SE. Estimates pooled across 20 imputed datasets. τ00 = random intercept variance. σ² = residual variance. ICC = intraclass correlation coefficient. Bold values indicate statistical significance (p < .05).
Table 3
Pooled Multilevel Model Estimates for Hostile/Rejection Parenting
Predictor
Model 1
Model 2
Model 3
Model 4
β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p β [95% CI] SE p
Fixed Effects
Intercept 0.126 [-0.065, 0.317] 0.097 0.195 0.135 [-0.057, 0.326] 0.098 0.169 0.135 [-0.057, 0.327] 0.098 0.168 0.145 [-0.045, 0.335] 0.097 0.136
SDQ Total 0.425 [0.285, 0.565] 0.071 < .001 0.423 [0.280, 0.566] 0.073 < .001 0.423 [0.280, 0.565] 0.073 < .001 0.427 [0.285, 0.569] 0.073 < .001
Child Age -0.057 [-0.198, 0.083] 0.072 0.424 -0.060 [-0.201, 0.080] 0.072 0.400 -0.060 [-0.201, 0.080] 0.072 0.401 -0.075 [-0.215, 0.066] 0.072 0.299
Child Sex -0.269 [-0.548, 0.010] 0.143 0.059 -0.287 [-0.570, -0.004] 0.144 0.047 -0.287 [-0.570, -0.004] 0.144 0.047 -0.309 [-0.592, -0.026] 0.144 0.032
Maternal Education -0.005 [-0.147, 0.137] 0.072 0.943 -0.017 [-0.163, 0.129] 0.074 0.816 -0.020 [-0.171, 0.131] 0.077 0.796 -0.022 [-0.175, 0.131] 0.078 0.778
Maternal ACEs — — — 0.005 [-0.149, 0.159] 0.079 0.948 0.003 [-0.153, 0.160] 0.080 0.966 -0.005 [-0.161, 0.150] 0.079 0.947
Informal Settlement Index — — — 0.012 [-0.135, 0.160] 0.075 0.869 0.012 [-0.136, 0.159] 0.075 0.874 0.015 [-0.132, 0.162] 0.075 0.840
DSM Cross-Cutting Symptom — — — 0.079 [-0.075, 0.234] 0.079 0.315 0.078 [-0.079, 0.234] 0.080 0.332 0.076 [-0.080, 0.231] 0.079 0.340
Mother Social Connections — — — — — — -0.011 [-0.163, 0.142] 0.078 0.892 -0.029 [-0.183, 0.125] 0.079 0.711
Community Participation — — — — — — — — — 0.107 [-0.035, 0.249] 0.073 0.140
Perceived Mutual Aid — — — — — — — — — 0.029 [-0.123, 0.182] 0.078 0.704
Variance Components
τ00 (between-cluster) 0.006 — — 0.005 — — 0.005 — — 0.001 — —
σ² (within-cluster) 0.792 — — 0.786 — — 0.786 — — 0.778 — —
ICC 0.007 — — 0.006 — — 0.006 — — 0.001 — —
Note. β = Standardized coefficient. 95% CI computed as β ± 1.96 × SE. Estimates pooled across 20 imputed datasets. τ00 = random intercept variance. σ² = residual variance. ICC = intraclass correlation coefficient. Bold values indicate statistical significance (p < .05).

I will run residual diagnostics to assess model assumptions of normality and homoscedasticity.

In addition to the use of eigenvalues outlined in the previous step, I will compute VIFs on the fitted models to assess multicollinearity.

3.3 Sensitivity Analyses

I will use 500 iterations of Monte Carlo simulations to randomly select one youth per family. This will result in a set of 500 subsamples with independent observations.

Code 
#|echo: false

#### Running models with simulations for random resampling ####
          
          
        ### first creating a version of merged_focal with all variables standardized
        merged_focal_std <- merged_focal %>%
    mutate(across(c(
    crpbi_comunicativo, 
    crpbi_controlador, 
    crpbi_hostil_rechazo,
    mother_aces, 
    informal_settlement_index, 
    mother_dsmtotal,
    mother_sc, 
    commparticipation, 
    mutual_aid,
    adult_education,
    child_demo_age,
    sdq_total
  ), ~ as.numeric(scale(.))))

          set.seed(123)
          n_simulations <- 500
          
          sim_results <- vector("list", n_simulations)
          
          for(i in 1:n_simulations) {
            
            # Randomly select one child per family
            subsample <- merged_focal_std %>%
              group_by(record_id) %>%
              slice_sample(n = 1) %>%
              ungroup()
            
            # OLS parallel to MLM model 4 for each outcome
            sim_results[[i]] <- list(
              com  = coef(lm(crpbi_comunicativo  ~ sdq_total + child_demo_age + 
                               child_demo_gender + adult_education + mother_aces +
                               informal_settlement_index + mother_dsmtotal + 
                               mother_sc + commparticipation + mutual_aid, 
                             data = subsample)),
              ctrl = coef(lm(crpbi_controlador   ~ sdq_total + child_demo_age + 
                               child_demo_gender + adult_education + mother_aces +
                               informal_settlement_index + mother_dsmtotal + 
                               mother_sc + commparticipation + mutual_aid, 
                             data = subsample)),
              hos  = coef(lm(crpbi_hostil_rechazo ~ sdq_total + child_demo_age + 
                               child_demo_gender + adult_education + mother_aces +
                               informal_settlement_index + mother_dsmtotal + 
                               mother_sc + commparticipation + mutual_aid, 
                             data = subsample))
            )
          }
          
 # Extract coefficients for each outcome into a dataframe
          com_coefs  <- do.call(rbind, lapply(sim_results, function(x) x$com))  
          ctrl_coefs <- do.call(rbind, lapply(sim_results, function(x) x$ctrl))
          hos_coefs  <- do.call(rbind, lapply(sim_results, function(x) x$hos))
          
          # Convert to long format for plotting
          com_long  <- as.data.frame(com_coefs)  %>% pivot_longer(everything(), names_to = "predictor", values_to = "estimate") %>% mutate(outcome = "Warm/Communicative")
          ctrl_long <- as.data.frame(ctrl_coefs) %>% pivot_longer(everything(), names_to = "predictor", values_to = "estimate") %>% mutate(outcome = "Controlling")
          hos_long  <- as.data.frame(hos_coefs)  %>% pivot_longer(everything(), names_to = "predictor", values_to = "estimate") %>% mutate(outcome = "Hostile/Rejecting")
          
          all_coefs <- bind_rows(com_long, ctrl_long, hos_long)

I estimated OLS regression models parallel to the ones outlined in the previous step with all the subsamples. I examined the distribution of the resulting fixed effect estimates for each of the outcomes across iterations, comparing their direction and magnitude in relation to an assumed null effect (0). These are presented in Figure 5.

I then examined the stability of the fixed effect estimates by comparing MLM coefficients from the full models against OLS coefficients derived from the 500 Monte Carlo subsamples. As can be seen in Table 1 MLM and OLS estimates for all predictors across all three parenting outcomes were close in proximity, suggesting that the results were not meaningfully influenced by the multilevel model structure or clustering assumptions.

Figure 5: Distribution of Fixed Effect Estimates across Iterations

Table 1: Sensitivity Analysis: OLS Monte Carlo Subsampling vs. MLM Pooled Estimates
Predictor
OLS Subsampling (n = 500)
MLM Pooled Estimate
Mean SD β SE p
Warm/Communicative
Child Age -0.093 0.062 -0.148 0.082 0.071
Child Sex -0.004 0.121 -0.106 0.163 0.518
Community Participation 0.022 0.053 0.033 0.083 0.688
DSM Cross-Cutting Symptom -0.094 0.085 -0.079 0.091 0.382
Informal Settlement Index -0.024 0.064 -0.060 0.086 0.489
Maternal ACEs 0.152 0.068 0.079 0.090 0.383
Maternal Education 0.055 0.052 0.051 0.088 0.564
Mother Social Connections 0.079 0.072 0.032 0.090 0.725
Perceived Mutual Aid -0.013 0.043 -0.049 0.088 0.581
SDQ Total -0.175 0.078 -0.187 0.084 0.026
Controlling
Child Age -0.298 0.039 -0.211 0.069 0.002
Child Sex -0.770 0.091 -0.559 0.133 < .001
Community Participation 0.070 0.046 0.052 0.085 0.539
DSM Cross-Cutting Symptom 0.099 0.054 0.074 0.096 0.439
Informal Settlement Index -0.008 0.047 0.037 0.088 0.671
Maternal ACEs 0.106 0.044 0.111 0.094 0.235
Maternal Education -0.080 0.044 -0.101 0.091 0.265
Mother Social Connections 0.021 0.039 0.028 0.092 0.759
Perceived Mutual Aid 0.007 0.035 -0.029 0.090 0.747
SDQ Total 0.281 0.045 0.324 0.072 < .001
Hostile/Rejecting
Child Age -0.107 0.053 -0.075 0.074 0.316
Child Sex -0.330 0.095 -0.308 0.149 0.038
Community Participation 0.126 0.036 0.109 0.077 0.157
DSM Cross-Cutting Symptom 0.130 0.076 0.078 0.084 0.352
Informal Settlement Index 0.024 0.059 0.015 0.079 0.849
Maternal ACEs -0.044 0.061 -0.007 0.084 0.934
Maternal Education -0.032 0.054 -0.021 0.082 0.795
Mother Social Connections -0.018 0.059 -0.028 0.083 0.736
Perceived Mutual Aid 0.033 0.039 0.029 0.082 0.724
SDQ Total 0.404 0.060 0.424 0.076 < .001
Note. OLS coefficients estimated across 500 Monte Carlo subsamples, each randomly selecting one child per family. β = standardized MLM coefficient from the full model refitted using Restricted Maximum Likelihood (REML), pooled across 20 imputed datasets. Bold values indicate statistical significance (p < .05).