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 (high warmth, low hostility) for positive youth mental health 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
While examining missingness per item, 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 160 youth nested within 100 mothers.
There were minimal data missing at the item-level (< 5% on all items). Little’s MCAR test indicated data were missing completely at random (MCAR) for youth reports (p = .181) but not for mother reports (p < .001). However, further examination of missing data patterns suggested missingness was clustered among a small number of cases and with no clear pattern. Given the results of Little’s MCAR test and the lack of a consistent pattern upon further examination, I assumed the data were likely missing at random (MAR) and proceeded to compute composite-level variables. Prior to computing the scales, I had determined that scale scores would be treated as missing if more than 20% of constituent items were missing.
Missingness at the composite-level was minimal, with 94.4% cases complete, and occurred predominantly at Level 2 (mother-level variables), ranging from 0.6% to 3.8% across variables. Table 1 presents a detailed view of missing values at the composite level.
Table 1: Missing Values at the Composite-Level
Variable
N Missing
% Missing
Biweekly Income
6
3.8
Warm/Communicative Parenting
2
1.2
Maternal Education
2
1.2
Hostile/Rejecting Parenting
1
0.6
Maternal ACEs
1
0.6
Informal Settlement Index
1
0.6
Mother DSM 5 Total
1
0.6
Mother Social Connections
1
0.6
Perceived Mutual Aid
1
0.6
Perceived Community Participation
1
0.6
1.1.3 Descriptives
After computing all composite variables, I examined measures of central tendency and dispersion for all observed focal variables. A summary of these can be found in Table 2. I also relied on boxplots (Figure 1, Figure 2) to examine the distribution of the data and screen for the possible presence of outliers. Through visual inspection of these plots, there was one youth whose age was outside of the target for this sample (19 years-old) but I decided to keep this observation given the proximity to the target age range and the absence of extreme values on key study variables. Beyond this case, I did not think any observations represented meaningful outliers for any of the study variables. There are a few aspects that are worth explictly noting. While controlling parenting was evenly spread, warm/communicative and hostile/rejecting parenting seemed to present with potential ceiling and floor effects, respectively. This limited range could attenuate associations otherwise present in the population.
At the family level, most variables showed relatively even speards. However, maternal ACEs and mother mental health symptoms showed positive skew whereas social connections showed a negative skew. Given that positive skew on mental health indicators is common in community samples, that the skew did not appear extreme, and that the multilevel estimator is generally considered robust to moderate non-normality, I retained the original metrics for all variables. As in with the parenting variables,it is worth noting that the restricted variance in these variables may attenuate associations.
Table 2: 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 externalizing
160
7.05
3.52
0.00
16.00
SDQ internalizing
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.
Figure 1: Boxplots to Examine Level 1 Variable Distribution
Figure 2: Boxplots to Examine Level 2 Variable Distribution
1.2 Multiple Imputation
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. 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; see supplementary imputation output below).
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
fml <-list( crpbi_controlador + crpbi_comunicativo + crpbi_hostil_rechazo ~1+ sdq_externalizing + sdq_internalizing + 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_externalizing, record_id) sdq_int_gpm <-clusterMeans(sdq_internalizing, record_id)})# Getting means from each dataset and summaring across themmeans_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 = variablesmeans_df <-do.call(rbind, means_list)rownames(means_df) <-paste0("imp_", 1:20)# Summary across all 20 imputationsround(means_df, 3)# And the overall average across all 20round(colMeans(means_df), 3)### Checking youth variablesyouth_vars <-c("sdq_externalizing", "sdq_internalizing", "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 datasetsimputed_means <-colMeans(means_df)imputed_means_df <-data.frame(variable =names(imputed_means),imp_mean =as.numeric(imputed_means))# Pivot to longmeans_long <-as.data.frame(means_df) %>%mutate(imp =1:20) %>%pivot_longer(-imp, names_to ="variable", values_to ="mean")# Mergeplot_df <-left_join(means_long, imputed_means_df, by ="variable")# Label mappingvar_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_dfplot_df <- plot_df %>%mutate(variable = dplyr::recode(variable, !!!var_labels))# Update imputed means df tooimputed_means_df <- imputed_means_df %>%mutate(variable = dplyr::recode(variable, !!!var_labels))# Re-mergeplot_df <-left_join( plot_df %>%select(-imp_mean), imputed_means_df, by ="variable")plot_df <- plot_df %>%mutate(mean_centered = mean - imp_mean)
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 3: Examining the Stability of Variable Means Across 20 Imputed Datasets
2 RQ 1: Parenting Dimensions and Youth Internalizing and Externalizing Behaviors
The goal of this first research question was to evaluate the role of parenting dimensions as promotive factors for youth mental health within the context of informal settlements. Specifically, whether the association between youths’ reports of caregiving practices and internalizing and externalizing symptoms mimic those found in the literature.
The hypotheses I had established were:
Hypothesis 1: Higher levels of parental warmth will be associated with lower levels of youth internalizing and externalizing symptoms.
Hypothesis 2: Higher levels of hostile and rejecting parenting will be associated with higher levels of youth internalizing and externalizing symptoms.
Hypothesis 3 (Exploratory): Controlling parenting will be associated with youth internalizing and externalizing symptoms, with the direction of these associations to be determined empirically.
As a useful reminder, to assess parenting practices for this study, I relied on a Spanish abbreviated version of the Child’s Report of Parent Behavior Inventory (CRPBI). Table 3 presents each of the dimensions included in the study with its corresponding items.
Table 3: Parenting Dimensions and Items
Warmt/Communicative
Hostile/Rejecting
Controlling
• Does he/she like to talk with you and tell you things? • Does he/she like to do things with you at home? • Does he/she speak to you with a soft and kind voice? • Do you feel better after telling him/her your problems? • Does he/she understand you when you tell him/her your problems? • Does he/she listen to your ideas and what you think? • Do you go together to interesting places and talk about the things you see there?
• Does he/she dislike the way you do things at home? • Does he/she say words that make you feel bad, like “stupid” or “dumb”? • Does he/she get very angry with you when you don’t help at home? • Does he/she get mad when you make noise in the house? • Does he/she act as if you bother him/her?
• Does he/she set a lot of rules in the house to maintain order? • Does he/she repeatedly tell you how you should do your chores/homework? • Does he/she want to control everything you do? • Does he/she try to change you? • Does he/she remind you of the things that are not allowed?
To test these hypotheses, I specified two sets of multilevle models. First, I specified three models for each internalizing and externalizing symptoms regressed on the three parenting dimensions. After examining these bivariate associations, I included youth age and youth sex in these models. The formal model equations for the unconditional and these two set of models follow.
\(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)
2.1 Regression Assumptions
I conducted a series of analyses to identify any potential violations of the regression assumptions of linearity, normality,homoskedasticity, and multicollinearity.
2.1.1 Linear Relationship
I examined the linearity of the relationship between each parenting predictor and the outcome variables using scatterplots with LOESS curves and estimated regression lines see (Figure 4 and Figure 5). Generally, the relationships between parenting dimensions and SDQ scores appeared approximately linear, with the LOESS curves not demonstrating any egregious violations of a linear fit and the confidence envelope of the LOESS curves encompassing the regression line, supporting the use of linear models. One potential exception is worth noting: the relationship between controlling parenting and externalizing symptoms, which showed a mild curvilinear pattern in the LOESS curve. Given the wide confidence envelope in this region, which still encompassed the regression line, I decided to retain linear models.
Figure 4: Scatterplot of Relationship between Parenting Dimensions and Internalizing
Figure 5: Scatterplot of Relationship between Parenting Dimensions and Externalizing
2.1.2 Normality
I examined the normality of model level 1 residuals visually using Q-Q plots and histograms of the residuals, and formally using the Shapiro-Wilk test. Residuals were approximately normally distributed, with no severe departures observed. Shapiro-Wilk tests indicated that residuals were normally distributed across all models (W range: 0.983–0.995), with the exception of the model predicting externalizing symptoms from communicative parenting, which showed marginal violations (W = 0.983, p = .047). Given the sensitivity of the Shapiro-Wilk test to sample size and the absence of meaningful deviations in the corresponding Q-Q plots, I did not interpret this as a substantive violation of the normality assumption.
Figure 6: Q-Q Plots for Examining Normality of the Level 1 Residuals for RQ1
Figure 7: Histograms for Examining Normality of the Level 1 Residuals for RQ1
I additionally examined the normality of Level 2 random effects (empirical Bayes residuals for the random intercept) using Q-Q plots, histograms, and Shapiro-Wilk tests, following @cho2022’s recommendations for multilevel models.
Random intercepts were approximately normally distributed across internalizing models (W range: 0.988–0.992, ps > .05). There were mild departures in normality for externalizing models in the Shapiro-Wilk tests (W range: 0.968–0.976, ps = .015–.061), but my visual inspection of histograms revealed only slight skew without extreme deviations.
Given that the data only had 100 clusters available at Level 2, the Shapiro-Wilk test has limited reliability and @cho2022 recommends visual inspection as the more informative criterion with a small number of clusters. Because of these considerations, I did not interpret these as substantive violations.
I examined homoskedasticity at Level 1 only, as all models included only a random intercept at Level 2 (i.e., no Level 2 predictors. As can be seen in Figure 10 assumption that residuals have constant variance looked reasonably met across all 6 models, with residuals scattered fairly randomly around zero without any obvious funneling or systematic patterns.
Code
residual_spread_rq1
Figure 10: Examining the Distribution of Model Residuals to Assess Homoskedasticity
2.1.4 Multicollinearity
I examined bivariate correlations among the predictor variables, computed variance inflation factors (VIF) and summed the reciprocal of the eigenvalues for all predictors included in Models 2. As seen in the heatmap below, the correlations among predictors were low to moderate (all |r| < .40). VIF values were all below the conventional threshold of 5 (range: 1.01–1.16), as can be seen in the output below.
Finally, I examined the sum of the reciprocals of the eigenvalues of the matrix of predictors for this second research question, and the total was 6.28, which provided additional evidence in support of there being no collinearity issues.
Model 1 examined the bivariate association between communicative parenting and each outcome, controlling for family clustering. Warm/communicative parenting was not a significant predictor of internalizing symptoms in Model 1 (β = -0.079, SE = 0.078, p = .314), nor did it reach significance in Model 2 after adding child age and child sex (β = -0.092, SE = 0.075, p = .219).
For externalizing symptoms, communicative parenting was a significant negative predictor in Model 1 (β = -0.169, SE = 0.078, p = .031), such that youth reporting higher levels of warm/communicative parenting showed lower externalizing symptoms. This effect remained significant and stable in Model 2 (β = -0.171, SE = 0.079, p = .030).
In terms of child covariates, both youth age and sex emerged as significant negative predictors of internalizing symptoms (β = -0.229, SE = 0.074, p = .002; β = 0.382, SE = 0.145, p = .008). This indicated that older and male youth reported lower internalizing symptoms.
For externalizing symptoms neither youth covariate was significant. However, child sex showed a near-significant trend (β = -0.268, SE = 0.153, p = .079), suggesting the possibility that female youth may report lower externalizing symptoms in this sample.
These results provide partial support for Hypothesis 1 — warm/communicative parenting was associated with lower externalizing symptoms as predicted, but the association with internalizing symptoms was not significant.
Code
tbl_comunicativo
Table 4: Pooled Multilevel Model Estimates for Warm/Communicative Parenting as a Predictor of Internalizing and Externalizing Behavior
Model 1
Model 2
Model 1
Model 2
Predictor
Internalizing Symptoms
Externalizing Symptoms
β [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).
2.2.2 Controlling Parenting
Controlling parenting was a significant positive predictor of both internalizing (β = 0.325, SE = 0.075, p < .001) and externalizing symptoms (β = 0.182, SE = 0.078, p = .020), such that youth reporting higher levels of controlling parenting showed higher levels of both internalizing and externalizing symptoms.
In Model 2, after adding child age and child sex, the effect of controlling parenting strengthened slightly for internalizing (β = 0.384, SE = 0.076, p < .001) and remained stable for externalizing (β = 0.191, SE = 0.083, p = .021).
Child sex emerged as a significant predictor of internalizing symptoms in Model 2 (β = 0.601, SE = 0.139, p < .001), indicating that female youth reported higher internalizing symptoms when accounting for controlling parenting and age. Child age showed a near-significant negative trend for internalizing (β = -0.129, SE = 0.071, p = .069). These results provide support for Hypothesis 3 — controlling parenting was associated with both internalizing and externalizing symptoms. In this sample, the direction of the association was positive across both outcomes.
Code
tbl_controlador
Table 5: Pooled Multilevel Model Estimates for Controlling Parenting as a Predictor of Internalizing and Externalizing Behavior
Model 1
Model 2
Model 1
Model 2
Predictor
Internalizing Symptoms
Externalizing Symptoms
β [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).
2.2.3 Hostile/Rejecting Parenting
Code
tbl_rechazo
Table 6: Pooled Multilevel Model Estimates for Hostile/Rejecting Parenting as a Predictor of Internalizing and Externalizing Behavior
Model 1
Model 2
Model 1
Model 2
Predictor
Internalizing Symptoms
Externalizing Symptoms
β [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).
Hostile/rejecting parenting was a significant positive predictor of both internalizing (β = 0.248, SE = 0.076, p = .001) and externalizing symptoms (β = 0.393, SE = 0.073, p < .001), such that youth reporting higher levels of hostile and rejecting parenting showed higher levels of both internalizing and externalizing symptoms.
In Model 2, after adding child age and child sex, the effect of hostile/rejecting parenting remained stable and significant for both internalizing (β = 0.256, SE = 0.072, p < .001) and externalizing (β = 0.395, SE = 0.073, p < .001).
Child age emerged as a significant negative predictor of internalizing symptoms in Model 2 (β = -0.192, SE = 0.072, p = .007), with older youth reporting lower internalizing symptoms. Child sex was also significant for internalizing (β = 0.459, SE = 0.140, p = .001), indicating that female youth reported higher internalizing symptoms when accounting for hostile/rejecting parenting and age. Neither child age nor child sex significantly predicted externalizing symptoms in Model 2.
These results provide support for Hypothesis 2 — hostile and rejecting parenting was associated with higher internalizing and externalizing symptoms in this sample, as predicted.
2.2.4 Overall Findings
Across all three parenting dimensions, the pattern of results provided at least partial support for the hypotheses. Hypothesis 1, which predicted that higher levels of parental warmth would be associated with lower internalizing and externalizing symptoms, was partially supported — warm/communicative parenting significantly predicted lower externalizing symptoms but was not a significant predictor of internalizing symptoms.
Hypothesis 2, which predicted that hostile and rejecting parenting would be associated with higher internalizing and externalizing symptoms, was fully supported — hostile/rejecting parenting was a significant positive predictor of both outcomes across both models. Hypothesis 3, which proposed that controlling parenting would be associated with youth symptoms with the direction to be determined empirically, was also supported — controlling parenting was a significant predictor of both internalizing and externalizing symptoms across both models. Importantly, this was a positive association so that higher levels of controlling parenting predicted higher levels of both internalizing and externalizing symptoms. These effects were of higher magnitude once child age and sex were included in the models.
Notably, hostile/rejecting parenting had a stronger effect on youth externalizing in comparison to internalizing symptoms, while the effect of controlling parenting was of higher magnitude for youth internalizing symptoms.
Child demographic characteristics emerged as consistent covariates for internalizing but not externalizing symptoms. Specifically, female and younger youth were more likely to report higher internalizing symptoms, once parenting dimensions were accounted for. The ICC values across models ranged from 0.040 to 0.280, indicating meaningful family-level clustering in youth outcomes, particularly for internalizing symptoms.
Model comparisons using the D3 pooled likelihood ratio test (@meng1992) indicated that adding child age and sex in Model 2 significantly improved model fit for internalizing symptoms across all three parenting dimensions: warm/communicative parenting (F(2) = 8.23, p < .001), controlling parenting (F(2) = 11.15, p < .001), and hostile/rejecting parenting (F(2) = 8.88, p < .001). In contrast, the addition of child age and sex did not significantly improve model fit for externalizing symptoms across any parenting dimension: warm/communicative parenting (F(2) = 1.82, p = .162), controlling parenting (F(2) = 1.68, p = .187), and hostile/rejecting parenting (F(2) = 1.90, p = .150) — suggesting that child demographic characteristics accounted for unique variance in internalizing but not externalizing symptoms beyond parenting.
3 RQ2: Identifying Risk factors and Promotive Factors for Mothers’ Parental Resilience in Informal Settlements
I drew on @gavidiapayne2015 to define parental resilience as “the capacity of parents to deliver competent, quality parenting to children despite adverse personal, family, and social circumstances.” (p. 113). Specifically, parental resilience regarding positive youth mental health outcomes. Drawing on RQ1 results, this meant high warm/communicative parenting, and low controlling and hostile/rejecting parenting.
I examined individual, relational, and community-level factors that either hindered (risk factors) or supported (promotive factors) mothers’ parental resilience in this context.
Hypothesis 4: Individual-level factors will predict parental resilience (good parenting quality) after controlling for youth characteristics and family demographics. Specifically, lower maternal mental health symptoms, higher educational attainment, age, will be associated with higher parental resilience (higher warmth, lower hostility), with the direction for controlling parenting to be determined empirically.
Hypothesis 5: Relational-level factors will predict parental resilience above and beyond individual-level factors and covariates. Specifically, higher social connections will be associated with higher parental resilience (higher warmth, lower hostility), with the direction for controlling parenting to be determined empirically. This tests whether relational resources explain unique variance in parenting quality beyond individual characteristics and adversity exposure.
Hypothesis 6: Community-level factors will predict parental resilience above and beyond individual and relational factors. Specifically, higher perceived community mutual aid and/or participation will be associated with higher parental resilience (higher warmth, lower hostility), with the direction for controlling parenting to be determined empirically. This tests whether community-level resources contribute to parenting quality beyond individual and relational determinants.
Hypothesis 7: According to resilience frameworks that consider both risks and protective factors, every ecological level is expected to make a unique contribution to parenting quality. More specifically, relational factors (Model 3) and community factors (Model 4) should add significant predictive value for parental resilience beyond what individual factors alone (Model 2) can explain. This will be shown through likelihood ratio tests, and greater reductions in Level 2 variance.
Following these hypotheses, four sequential models tested unique variance explained by each ecological level. The goal with this approach was to examine whether more distal factors (e.g., 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.
All Level 2 predictors were grand mean centered and the models were estimated employing Maximum Likelihood Estimation (MLE) using the lme4 package in R (@bates2015) 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 (@finch2019). I used these final estimates to compare with the standard OLS estimates in the final step of analyses. Formal model equations follow.
3.1 Regression Assumptions
I conducted a series of analyses to identify any potential violations of the regression assumptions of linearity, normality, homoskedasticity, and multicollinearity.
3.1.1 Linear Relationship
I examined Level 1 linearity by plotting marginal standardized residuals against each continuous Level 1 predictor, following @cho2022) (see Figure 11 and Figure 12). Residuals were generally scattered around zero across the full range of both SDQ total problems and child age for all three parenting outcomes, with no systematic nonlinear patterns observed, supporting the assumption of linearity for Level 1 predictors. A small number of observations showed marginal standardized residuals beyond ±3, most notably in the Communicative Parenting model where a few observations approached ±4. Given the bounded 1–3 scale of the parenting outcomes, I attributed these extreme residuals to the discrete distributional constraints of the outcome rather than true model violations, and noted them as a limitation of applying linear models to bounded Likert-type outcomes. I assessed Level 2 linearity via plots of standardized EB residuals against each Level 2 predictor (see Figure 18, Figure 19, and Figure 20), which showed no systematic patterns inconsistent with linearity.
Figure 11: Level 1 Linearity: Marginal Standardized Residuals vs SDQ Total Problems (RQ2)
Figure 12: Level 1 Linearity: Marginal Standardized Residuals vs Child Age (RQ2)
3.1.2 Normality
I examined the normality of model Level 1 residuals visually using Q-Q plots and histograms of the residuals, and formally using the Shapiro-Wilk test. Shapiro-Wilk tests indicated significant departures from normality for the Warm/Communicative (W = 0.929, p < .001) and Hostile/Rejecting models (W = 0.974, p = .004), while the Controlling model showed no significant violation (W = 0.986, p = .103). Visual inspection of the histograms revealed left skew across all three models, most pronounced for Warm/Communicative parenting, which is consistent with the ceiling effect observed in the descriptive statistics for this variable. Given that these departures are attributable to the bounded 1–3 scale of the parenting outcomes rather than model misspecification, and that multilevel models are generally robust to moderate violations of normality with respect to fixed effect estimates (as discussed by @cho2022), I did not interpret these as substantive violations.
Figure 13: Q-Q Plots for Examining Normality of the Level 1 Residuals for RQ2
Figure 14: Histograms for Examining Normality of the Level 1 Residuals for RQ2
Model W p
1 Warm/Communicative M4 0.929 0.000
2 Controlling M4 0.986 0.103
3 Hostile/Rejecting M4 0.974 0.004
I additionally examined the normality of Level 2 random effects (empirical Bayes residuals for the random intercept) using Q-Q plots, histograms, and Shapiro-Wilk tests, following @cho2022’s recommendations for multilevel models. Shapiro-Wilk tests indicated that residuals were normally distributed across all models (W range: .977–.993), with no models showing significant violations (all ps > .05).
Figure 15: Q-Q Plots of Level 2 Empirical Bayes Residuals (Random Intercepts) for RQ2
Figure 16: Histograms of Level 2 Empirical Bayes Residuals (Random Intercepts) for RQ2
Model W p
1 Warm/Communicative M4 0.977 0.079
2 Controlling M4 0.992 0.843
3 Hostile/Rejecting M4 0.993 0.861
3.1.3 Homoskedasticity
I examined homoskedasticity at both Level 1 and Level 2, as models included Level 2 predictors. As can be seen in Figure 17, the assumption that residuals have constant variance looked reasonably met across all three models, with residuals scattered fairly randomly around zero without any obvious funneling or systematic patterns. At Level 2, I examined whether the variance of the standardized EB residuals was constant across the range of each Level 2 predictor (see Figure 18, Figure 19, and Figure 20). The Hostile/Rejecting model showed no substantive concerns. The Communicative Parenting model showed mild curvature for some predictors, most evident at the extremes of the distributions where data are sparse. The Controlling Parenting model showed the most pronounced patterns, particularly for perceived mutual aid and social connections. These patterns are noted as a limitation in the interpretation of findings for controlling parenting specifically.
Figure 17: Examining the Distribution of Model Residuals to Assess Level 1 Homoskedasticity for RQ2
Figure 18: Level 2 Heteroskedasticity: Standardized EB Residuals vs Level 2 Predictors for Communicative Parenting
Figure 19: Level 2 Heteroskedasticity: Standardized EB Residuals vs Level 2 Predictors for Controlling Parenting
Figure 20: Level 2 Heteroskedasticity: Standardized EB Residuals vs Level 2 Predictors for Hostile/Rejecting Parenting
3.1.4 Multicollinearity
I examined bivariate correlations among the predictor variables, computed variance inflation factors (VIF), and summed the reciprocal of the eigenvalues for all predictors included in Model 4. As seen in the heatmap below, the correlations among predictors were low to moderate (all |r| < .50). VIF values were all below the conventional threshold of 5 (range: 1.05–1.29), as can be seen in the output below.
Finally, I examined the sum of the reciprocals of the eigenvalues of the matrix of predictors for this research question, and the total was 12, which provided additional evidence in support of there being no collinearity issues.
Model 1 included individual-level predictors: youth total problems, child age, child sex, and maternal education. SDQ total problems significantly predicted communicative parenting (β = -0.174, SE = 0.079, p = .028), such that youth with higher levels of total problems perceived themselves as receiving less warm and communicative parenting. Child age also emerged as a significant predictor (β = -0.156, SE = 0.079, p = .048), with older youth reporting lower levels of communicative parenting. Child sex and maternal education were not significant predictors of communicative parenting.
Model 2 added maternal adversity indicators (maternal ACEs, informal settlement index, and DSM cross-cutting symptoms), none of which significantly predicted communicative parenting (all ps > .30). The effects of SDQ total (β = -0.184, SE = 0.081, p = .023) remained significant, while child age became marginally non-significant (β = -0.146, SE = 0.079, p = .063).
Model 3 added mother social connections, which was not a significant predictor (β = 0.024, SE = 0.084, p = .773). Model 4 added community-level predictors (community participation and perceived mutual aid), neither of which significantly predicted communicative parenting (all ps > .55). In the final model, only SDQ total problems remained a significant predictor (β = -0.188, SE = 0.081, p = .021), with child age remaining marginally non-significant (β = -0.149, SE = 0.079, p = .060).
None of the mother-level or community-level predictors reached significance across any model. These findings are summarized in Table 7.
Code
tbl_com
Table 7: Pooled Multilevel Model Estimates for Warm/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).
3.2.2 Controlling Parenting
Model 1 included youth total problems, child age, child sex, and maternal education as predictors. SDQ total problems significantly predicted controlling parenting (β = 0.324, SE = 0.069, p < .001), such that youth with higher total problems reported higher levels of controlling parenting. Child age (β = -0.213, SE = 0.067, p = .002) and child sex (β = -0.564, SE = 0.129, p < .001) were also significant, with older youth and female youth reporting lower levels of controlling parenting. Maternal education was not significant (β = -0.080, SE = 0.081, p = .321).
Model 2 added maternal adversity indicators, none of which significantly predicted controlling parenting (all ps > .23). The effects of SDQ total, child age, and child sex remained stable and significant across Models 2 and 3 with minimal change in magnitude. Model 3 added mother social connections, which was not a significant predictor (β = 0.027, SE = 0.086, p = .756). Model 4 added community-level predictors, neither of which significantly predicted controlling parenting (all ps > .50).
In the final model, SDQ total (β = 0.323, SE = 0.070, p < .001), child age (β = -0.214, SE = 0.067, p = .001), and child sex (β = -0.567, SE = 0.130, p < .001) remained significant predictors, with no mother-level or community-level predictors reaching significance across any model. These findings are summarized in Table 8.
Code
tbl_ctrl
Table 8: 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).
3.2.3 Hostile/Rejecting Parenting
Model 1 included youth total problems, child age, child sex, and maternal education as predictors. SDQ total problems was a significant predictor of hostile/rejecting parenting (β = 0.425, SE = 0.071, p < .001), such that youth with higher total problems reported higher levels of hostile and rejecting parenting. Child sex showed a marginal trend in Model 1 (β = -0.269, SE = 0.143, p = .059), suggesting female youth may report lower hostile/rejecting parenting. Child age and maternal education were not significant predictors.
Model 2 added maternal adversity indicators, none of which significantly predicted hostile/rejecting parenting (all ps > .30). Child sex reached significance in Model 2 (β = -0.287, SE = 0.144, p = .047) and remained significant in Model 3. SDQ total remained stable and significant across all models. Model 3 added mother social connections, which was not a significant predictor (β = -0.011, SE = 0.078, p = .892). Model 4 added community-level predictors, with community participation showing a near-significant trend in the opposite direction than hypothesized (β = 0.107, SE = 0.073, p = .140) and perceived mutual aid not significant (β = 0.029, SE = 0.078, p = .704).
In the final model, SDQ total (β = 0.427, SE = 0.073, p < .001) and child sex (β = -0.309, SE = 0.144, p = .032) were significant predictors, with no mother-level or community-level predictors reaching significance.
Code
tbl_hos
Table 9: Pooled Multilevel Model Estimates for Hostile/Rejecting 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).
3.2.4 Overall Findings
Across all three parenting outcomes none of the mother-level or community-level predictors reached statistical significance across any outcome or model, providing no support for Hypotheses 4 through 6 regarding the role of individual maternal characteristics, relational resources, or community factors in predicting parenting quality. These null findings should be interpreted with caution given the study was underpowered to detect small effects at Level 2, with only 100 mothers at the cluster level.
It is also worth noting the ICC pattern across outcomes — controlling parenting showed substantially higher between-family clustering (ICC = 0.355–0.382) compared to communicative and hostile/rejecting parenting (ICC < 0.01), suggesting that controlling parenting behavior is more consistent within families across siblings than the other parenting dimensions. In other words, warm/communicative and hostile/rejecting parenting seems to be more influenced by youth-level characteristics whereas controlling parenting seems to be driven by family-level factors, albeit none of those included in this study being significant. This might be a function of the underpowered nature of this study, measurement error, or the need to identify other family-level factors that could explain between family variance in controlling parenting.
Sequential model comparisons using the D3 pooled likelihood ratio test indicated that no block of predictors significantly improved model fit beyond the previous model for any parenting outcome. These results indicate that neither individual maternal characteristics, relational resources, and community-level factors explained significant additional variance in parenting quality beyond youth characteristics, again, providing no support for Hypotheses 4 through 7.
3.3 Sensitivity Analyses
I used 500 iterations of Monte Carlo simulations to randomly select one youth per family. This resulted 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 in1: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 21.
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 10 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, clustering assumptions, or normality violations observed at Level 1.
Given that the study was underpowered to detect small effects at Level 2 due to the sample size of 100 mothers, non-significant results should be viewed with caution. The OLS subsampling analyses corroborate the directional consistency of the estimates; however, the substantial standard deviations observed for several predictors highlight the expected instability in the estimates associated with a limited Level 2 sample.
Figure 21: Distribution of OLS Coefficients Across 500 Monte Carlo Subsamples
Table 10: 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).
Source Code
---title: "Study 2 Results"author: "Ana M. Díaz"format: html: embed-resources: true code-fold: true code-tools: true toc: true toc-float: trueexecute: echo: true results: false warning: falseeditor: sourcetoc-location: leftfig-cap-location: topnumber-sections: true---```{=html}<style> body { font-family: "Times New Roman", Times, serif; }</style>``````{r}#| label: prep#| echo: false#| include: false#### Loading libraries and data ####library(ggplot2)library(patchwork)library(dplyr)library(tidyverse)library(Hmisc)library(haven)library(patchwork)library(tidyr)library(gtsummary)library(stringr)library(RColorBrewer)library(reshape2)library(ggmosaic)library(kableExtra)library(sjlabelled)library(lme4)library(performance)library(MLMusingR)library(jomo)library(mitml)library(gt)library(corrplot)library(car)library(polycor)library(HLMdiag)load("~/Dropbox/Research/Informal settlements/Datos y Analisis/analisis-encuestas/data/asentamientos.2025.2.15.RData")```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 (high warmth, low hostility) for positive youth mental health 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?## Preliminary Steps### Computing and Examining Focal Variables#### Item-Level MissingnessWhile examining missingness per item, 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 160 youth nested within 100 mothers.There were minimal data missing at the item-level (\< 5% on all items). Little's MCAR test indicated data were missing completely at random (MCAR) for youth reports (p = .181) but not for mother reports (p \< .001). However, further examination of missing data patterns suggested missingness was clustered among a small number of cases and with no clear pattern. Given the results of Little's MCAR test and the lack of a consistent pattern upon further examination, I assumed the data were likely missing at random (MAR) and proceeded to compute composite-level variables. Prior to computing the scales, I had determined that scale scores would be treated as missing if more than 20% of constituent items were missing.```{r}#| label: prepping-the-data#| echo: true#| results: false#### First creating a function to remove youth data from mom data and viceversadrop_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 columnsrecode_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) elseNA_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 scoressdq_externalizing = sdq_conduct + sdq_hyperactivity,sdq_internalizing = 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 scoresace_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(###### Infrastructurewater_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 vulnerabilitiesflooding =ifelse(inundaciones >0, 1, 0),###### Tenure insecurityno_land_title =ifelse(terrenos2 ==0, 1, 0),eviction_exp =ifelse(terrenos4 ==1, 1, 0),###### Community safetycommunity_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 datarowwise() %>%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_externalizing, sdq_internalizing, 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 childmerged_focal <- youth_data_focal %>%left_join(mother_data_focal, by ="record_id")# Check structure looks rightmerged_focal %>%count(record_id) %>%summary() ```#### Composite-Level MissingnessMissingness at the composite-level was minimal, with 94.4% cases complete, and occurred predominantly at Level 2 (mother-level variables), ranging from 0.6% to 3.8% across variables. @tbl-missing presents a detailed view of missing values at the composite level.```{r}#| label: tbl-missing#| echo: false#| tbl-cap: "Missing Values at the Composite-Level"naniar::miss_var_summary(merged_focal) %>%filter(pct_miss >0) %>%mutate(pct_miss =round(pct_miss, 1)) %>%mutate(variable = dplyr::recode(variable,"quincena"="Biweekly Income","crpbi_comunicativo"="Warm/Communicative Parenting","adult_education"="Maternal Education","crpbi_hostil_rechazo"="Hostile/Rejecting Parenting","mother_aces"="Maternal ACEs","informal_settlement_index"="Informal Settlement Index","mother_sc"="Mother Social Connections","mutual_aid"="Perceived Mutual Aid","mother_dsmtotal"="Mother DSM 5 Total","commparticipation"="Perceived Community Participation" )) %>%rename(Variable = variable, "N Missing"= n_miss, "% Missing"= pct_miss) %>% knitr::kable(align =c("l", "c", "c"),booktabs =TRUE, format ="html" )```#### DescriptivesAfter computing all composite variables, I examined measures of central tendency and dispersion for all observed focal variables. A summary of these can be found in @tbl-descriptives. I also relied on boxplots (@fig-boxplot1, @fig-boxplot2) to examine the distribution of the data and screen for the possible presence of outliers. Through visual inspection of these plots, there was one youth whose age was outside of the target for this sample (19 years-old) but I decided to keep this observation given the proximity to the target age range and the absence of extreme values on key study variables. Beyond this case, I did not think any observations represented meaningful outliers for any of the study variables. There are a few aspects that are worth explictly noting. While controlling parenting was evenly spread, warm/communicative and hostile/rejecting parenting seemed to present with potential ceiling and floor effects, respectively. This limited range could attenuate associations otherwise present in the population.At the family level, most variables showed relatively even speards. However, maternal ACEs and mother mental health symptoms showed positive skew whereas social connections showed a negative skew. Given that positive skew on mental health indicators is common in community samples, that the skew did not appear extreme, and that the multilevel estimator is generally considered robust to moderate non-normality, I retained the original metrics for all variables. As in with the parenting variables,it is worth noting that the restricted variance in these variables may attenuate associations.```{r}#| echo: false#| warning: false#### descriptives ##### For continuous variablescontinuous_vars <-c("sdq_externalizing", "sdq_internalizing", "sdq_total","crpbi_comunicativo", "crpbi_controlador", "crpbi_hostil_rechazo","child_demo_age","mother_aces", "informal_settlement_index","mother_dsmtotal", "mother_sc", "mutual_aid", "commparticipation", "adult_education", "quincena")var_labels <-c(sdq_externalizing ="SDQ externalizing",sdq_internalizing ="SDQ internalizing",sdq_total ="SDQ Total",crpbi_comunicativo ="Communicative Parenting",crpbi_controlador ="Controlling Parenting",crpbi_hostil_rechazo ="Hostile/Rejecting Parenting",child_demo_age ="Child Age",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",adult_education ="Maternal Education",quincena ="Biweekly Income")group_labels <-c(sdq_externalizing ="Child Outcomes",sdq_internalizing ="Child Outcomes",sdq_total ="Child Outcomes",crpbi_comunicativo ="Parenting",crpbi_controlador ="Parenting",crpbi_hostil_rechazo ="Parenting",child_demo_age ="Child Demographics",mother_aces ="Mother Characteristics",informal_settlement_index ="Mother Characteristics",mother_dsmtotal ="Mother Characteristics",mother_sc ="Mother Characteristics",mutual_aid ="Community",commparticipation ="Community",adult_education ="Mother Characteristics",quincena ="Mother Characteristics")desc_l1 <- merged_focal %>%select(all_of(c("sdq_externalizing", "sdq_internalizing", "sdq_total","crpbi_comunicativo", "crpbi_controlador", "crpbi_hostil_rechazo","child_demo_age"))) %>%summarise(across(everything(), list(n =~sum(!is.na(.)),mean =~mean(., na.rm =TRUE),sd =~sd(., na.rm =TRUE),min =~min(., na.rm =TRUE),max =~max(., na.rm =TRUE) ))) %>%pivot_longer(everything(), names_to =c("variable", ".value"), names_sep ="_(?=[^_]+$)")desc_l2 <- merged_focal %>%distinct(record_id, .keep_all =TRUE) %>%select(all_of(c("mother_aces", "informal_settlement_index", "mother_dsmtotal","mother_sc", "mutual_aid", "commparticipation","adult_education", "quincena"))) %>%summarise(across(everything(), list(n =~sum(!is.na(.)),mean =~mean(., na.rm =TRUE),sd =~sd(., na.rm =TRUE),min =~min(., na.rm =TRUE),max =~max(., na.rm =TRUE) ))) %>%pivot_longer(everything(), names_to =c("variable", ".value"), names_sep ="_(?=[^_]+$)")desc_all <-bind_rows(desc_l1, desc_l2) %>%mutate(label = dplyr::recode(variable, !!!var_labels),group = dplyr::recode(variable, !!!group_labels),group =factor(group, levels =c("Child Demographics", "Child Outcomes", "Parenting", "Mother Characteristics", "Community")) ) %>%arrange(group, variable) %>%mutate(across(c(mean, sd, min, max), ~sprintf("%.2f", .)))### table ###descriptives_table <- desc_all %>%select(group, label, n, mean, sd, min, max) %>%gt(rowname_col ="label", groupname_col ="group") %>%cols_label(n ="n",mean ="M",sd ="SD",min ="Min",max ="Max" ) %>%tab_stubhead(label ="Variable") %>%tab_source_note(source_note ="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." ) %>%tab_options(table.border.top.style ="solid",table.border.top.width =px(2),table.border.bottom.style ="solid",table.border.bottom.width =px(2),heading.border.bottom.style ="solid",heading.border.bottom.width =px(1),column_labels.border.top.style ="none",column_labels.border.bottom.style ="solid",column_labels.border.bottom.width =px(1),table_body.border.top.style ="none",table_body.border.bottom.style ="none",stub.border.style ="none",row_group.border.top.style ="solid",row_group.border.top.width =px(1),row_group.border.bottom.style ="none",row.striping.include_table_body =FALSE,table.font.size =px(12),table.font.names ="Times New Roman",heading.title.font.size =px(12),heading.subtitle.font.size =px(12),heading.align ="left",data_row.padding =px(4) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_title("title") ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_title("subtitle") ) %>%tab_style(style =cell_text(align ="left"),locations =cells_stub() ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_row_groups() )``````{r}#| label: tbl-descriptives#| tbl-cap: "Descriptive Statistics for Observed Study Variables"#| echo: falsedescriptives_table``````{r}#| label: fig-boxplot1#| fig-cap: "Boxplots to Examine Level 1 Variable Distribution"#| echo: false youth_data_focal %>%select(-record_id, -child_demo_gender) %>%pivot_longer(everything(), names_to ="variable", values_to ="value") %>%mutate(variable = dplyr::recode(variable,"sdq_externalizing"="SDQ externalizing","sdq_internalizing"="SDQ internalizing","sdq_total"="SDQ Total","crpbi_comunicativo"="Communicative Parenting","crpbi_controlador"="Controlling Parenting","crpbi_hostil_rechazo"="Hostile/Rejecting Parenting","child_demo_age"="Child Age" )) %>%ggplot(aes(x = variable, y = value)) +geom_boxplot(outlier.shape =NA, fill ="lightgreen", alpha =0.5) +geom_jitter(width =0.15, alpha =0.3, size =0.8) +facet_wrap(~ variable, scales ="free") +theme_minimal() +theme(axis.text.x =element_blank(),axis.ticks.x =element_blank()) +labs(x =NULL, y =NULL) ``````{r}#| label: fig-boxplot2#| fig-cap: "Boxplots to Examine Level 2 Variable Distribution"#| echo: falsemother_data_focal %>%select(-record_id) %>%pivot_longer(everything(), names_to ="variable", values_to ="value") %>%mutate(variable = dplyr::recode(variable,"quincena"="Biweekly Income","mother_aces"="Maternal ACEs","informal_settlement_index"="Informal Settlement Index","mother_dsmtotal"="DSM Cross-Cutting","mother_sc"="Social Connections","mutual_aid"="Perceived Mutual Aid","commparticipation"="Community Participation","adult_education"="Maternal Education" )) %>%ggplot(aes(x = variable, y = value)) +geom_boxplot(outlier.shape =NA, fill ="lightblue", alpha =0.5) +geom_jitter(width =0.15, alpha =0.3, size =0.8) +facet_wrap(~ variable, scales ="free") +theme_minimal() +theme(axis.text.x =element_blank(),axis.ticks.x =element_blank()) +labs(x =NULL, y =NULL)```### Multiple ImputationI 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. 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; see supplementary imputation output below).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.```{r}#| label: multiple-imputation#| echo: true#| results: falsefml <-list( crpbi_controlador + crpbi_comunicativo + crpbi_hostil_rechazo ~1+ sdq_externalizing + sdq_internalizing + 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_externalizing, record_id) sdq_int_gpm <-clusterMeans(sdq_internalizing, record_id)})# Getting means from each dataset and summaring across themmeans_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 = variablesmeans_df <-do.call(rbind, means_list)rownames(means_df) <-paste0("imp_", 1:20)# Summary across all 20 imputationsround(means_df, 3)# And the overall average across all 20round(colMeans(means_df), 3)### Checking youth variablesyouth_vars <-c("sdq_externalizing", "sdq_internalizing", "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)``````{r}#| label: imputed-means#| echo: true#| results: false# Overall mean across the 20 imputed datasetsimputed_means <-colMeans(means_df)imputed_means_df <-data.frame(variable =names(imputed_means),imp_mean =as.numeric(imputed_means))# Pivot to longmeans_long <-as.data.frame(means_df) %>%mutate(imp =1:20) %>%pivot_longer(-imp, names_to ="variable", values_to ="mean")# Mergeplot_df <-left_join(means_long, imputed_means_df, by ="variable")# Label mappingvar_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_dfplot_df <- plot_df %>%mutate(variable = dplyr::recode(variable, !!!var_labels))# Update imputed means df tooimputed_means_df <- imputed_means_df %>%mutate(variable = dplyr::recode(variable, !!!var_labels))# Re-mergeplot_df <-left_join( plot_df %>%select(-imp_mean), imputed_means_df, by ="variable")plot_df <- plot_df %>%mutate(mean_centered = mean - imp_mean)```::: {.callout-note collapse="true"}## Imputation Model Output```{r}#| echo: falsesummary(imp)```:::```{r}#| label: fig-imputations#| fig-cap: "Examining the Stability of Variable Means Across 20 Imputed Datasets"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()``````{r}#| label: sdq-total-imputed#| echo: false#| results: false## computing SDQ total for imputed dataimplist <-within(implist, { sdq_total <- sdq_internalizing + sdq_externalizing})youth_data$sdq_total```## RQ 1: Parenting Dimensions and Youth Internalizing and Externalizing BehaviorsThe goal of this first research question was to evaluate the role of parenting dimensions as promotive factors for youth mental health within the context of informal settlements. Specifically, whether the association between youths’ reports of caregiving practices and internalizing and externalizing symptoms mimic those found in the literature.The hypotheses I had established were:**Hypothesis 1:** Higher levels of parental warmth will be associated with lower levels of youth internalizing and externalizing symptoms.**Hypothesis 2:** Higher levels of hostile and rejecting parenting will be associated with higher levels of youth internalizing and externalizing symptoms.**Hypothesis 3** (Exploratory)**:** Controlling parenting will be associated with youth internalizing and externalizing symptoms, with the direction of these associations to be determined empirically.As a useful reminder, to assess parenting practices for this study, I relied on a Spanish abbreviated version of the Child's Report of Parent Behavior Inventory (CRPBI). @tbl-parenting-items presents each of the dimensions included in the study with its corresponding items.| Warmt/Communicative | Hostile/Rejecting | Controlling ||---------------------|----------------------------|-----------------------|| • Does he/she like to talk with you and tell you things? <br> • Does he/she like to do things with you at home? <br> • Does he/she speak to you with a soft and kind voice? <br> • Do you feel better after telling him/her your problems? <br> • Does he/she understand you when you tell him/her your problems? <br> • Does he/she listen to your ideas and what you think? <br> • Do you go together to interesting places and talk about the things you see there? | • Does he/she dislike the way you do things at home? <br> • Does he/she say words that make you feel bad, like "stupid" or "dumb"? <br> • Does he/she get very angry with you when you don't help at home? <br> • Does he/she get mad when you make noise in the house? <br> • Does he/she act as if you bother him/her? | • Does he/she set a lot of rules in the house to maintain order? <br> • Does he/she repeatedly tell you how you should do your chores/homework? <br> • Does he/she want to control everything you do? <br> • Does he/she try to change you? <br> • Does he/she remind you of the things that are not allowed? |: Parenting Dimensions and Items {#tbl-parenting-items}To test these hypotheses, I specified two sets of multilevle models. First, I specified three models for each internalizing and externalizing symptoms regressed on the three parenting dimensions. After examining these bivariate associations, I included youth age and youth sex in these models. The formal model equations for the unconditional and these two set of models follow.::: {style="text-align: center"}**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*:::::: {style="margin-left: 2em"}**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):::### Regression AssumptionsI conducted a series of analyses to identify any potential violations of the regression assumptions of linearity, normality,homoskedasticity, and multicollinearity.#### Linear RelationshipI examined the linearity of the relationship between each parenting predictor and the outcome variables using scatterplots with LOESS curves and estimated regression lines see (@fig-rq-diagnostics-int and @fig-rq-diagnostics-ext). Generally, the relationships between parenting dimensions and SDQ scores appeared approximately linear, with the LOESS curves not demonstrating any egregious violations of a linear fit and the confidence envelope of the LOESS curves encompassing the regression line, supporting the use of linear models. One potential exception is worth noting: the relationship between controlling parenting and externalizing symptoms, which showed a mild curvilinear pattern in the LOESS curve. Given the wide confidence envelope in this region, which still encompassed the regression line, I decided to retain linear models.```{r}#| label: rq1-diagnostics#| echo: false### RUNNING RQ 1 REGRESSION DIAGNOSTICS ###### internalizing scatterplotsp1 <-ggplot(merged_focal, aes(x = crpbi_comunicativo, y = sdq_internalizing)) +geom_point(alpha =0.5) +geom_smooth(method ="loess", color ="blue") +geom_smooth(method ="lm", color ="red", linetype ="dashed") +labs(title ="Warm/Communicative & internalizing", x ="Warm/Communicative", y ="internalizing") +theme_minimal()p2 <-ggplot(merged_focal, aes(x = crpbi_controlador, y = sdq_internalizing)) +geom_point(alpha =0.5) +geom_smooth(method ="loess", color ="blue") +geom_smooth(method ="lm", color ="red", linetype ="dashed") +labs(title ="Controlling & internalizing", x ="Controlling", y ="internalizing")+theme_minimal()p3 <-ggplot(merged_focal, aes(x = crpbi_hostil_rechazo, y = sdq_internalizing)) +geom_point(alpha =0.5) +geom_smooth(method ="loess", color ="blue") +geom_smooth(method ="lm", color ="red", linetype ="dashed") +labs(title ="Hostile/Rejecting & internalizing", x ="Hostile/Rejecting", y ="internalizing")+theme_minimal()# externalizing scatterplotsp4 <-ggplot(merged_focal, aes(x = crpbi_comunicativo, y = sdq_externalizing)) +geom_point(alpha =0.5) +geom_smooth(method ="loess", color ="blue") +geom_smooth(method ="lm", color ="red", linetype ="dashed") +labs(title ="Warm/Communicative & externalizing", x ="Warm/Communicative", y ="externalizing")+theme_minimal()p5 <-ggplot(merged_focal, aes(x = crpbi_controlador, y = sdq_externalizing)) +geom_point(alpha =0.5) +geom_smooth(method ="loess", color ="blue") +geom_smooth(method ="lm", color ="red", linetype ="dashed") +labs(title ="Controlling & externalizing", x ="Controlling", y ="externalizing")+theme_minimal()p6 <-ggplot(merged_focal, aes(x = crpbi_hostil_rechazo, y = sdq_externalizing)) +geom_point(alpha =0.5) +geom_smooth(method ="loess", color ="blue") +geom_smooth(method ="lm", color ="red", linetype ="dashed") +labs(title ="Hostile/Rejecting & externalizing", x ="Hostile/Rejecting", y ="externalizing")+theme_minimal()implist_std <-within(implist, { sdq_internalizing <-scale(sdq_internalizing) sdq_externalizing <-scale(sdq_externalizing) sdq_total <-scale(sdq_total) crpbi_comunicativo <-scale(crpbi_comunicativo) crpbi_controlador <-scale(crpbi_controlador) crpbi_hostil_rechazo <-scale(crpbi_hostil_rechazo) child_demo_age <-scale(child_demo_age) mother_aces <-scale(mother_aces) informal_settlement_index <-scale(informal_settlement_index) mother_dsmtotal <-scale(mother_dsmtotal) mother_sc <-scale(mother_sc) commparticipation <-scale(commparticipation) mutual_aid <-scale(mutual_aid) adult_education <-scale(adult_education) quincena <-scale(quincena)})``````{r}#| label: fig-rq-diagnostics-int#| fig-cap: Scatterplot of Relationship between Parenting Dimensions and Internalizing#| fig-cap-location: top#| echo: false# Plot all together(p1 / p2 / p3)``````{r}#| label: fig-rq-diagnostics-ext#| fig-cap: "Scatterplot of Relationship between Parenting Dimensions and Externalizing"#| fig-cap-location: top#| echo: false# Plot all together(p4 / p5 / p6)``````{r}#| label: regression-diagnostics-rq1#| echo: true#| include: false#| results: falseone_imp <- implist_std[[1]] # i am using the 1/20 imputed datasets### fitting all final models with just this one datasetmodels_rq1 <-list("Int-Communicative + Covs"=lmer(sdq_internalizing ~ crpbi_comunicativo + child_demo_age + child_demo_gender + (1| record_id), data = one_imp, REML =FALSE),"Int-Controlling + Covs"=lmer(sdq_internalizing ~ crpbi_controlador + child_demo_age + child_demo_gender + (1| record_id), data = one_imp, REML =FALSE),"Int-Hostile + Covs"=lmer(sdq_internalizing ~ crpbi_hostil_rechazo + child_demo_age + child_demo_gender + (1| record_id), data = one_imp, REML =FALSE),"Ext-Communicative + Covs"=lmer(sdq_externalizing ~ crpbi_comunicativo + child_demo_age + child_demo_gender + (1| record_id), data = one_imp, REML =FALSE),"Ext-Controlling + Covs"=lmer(sdq_externalizing ~ crpbi_controlador + child_demo_age + child_demo_gender + (1| record_id), data = one_imp, REML =FALSE),"Ext-Hostile + Covs"=lmer(sdq_externalizing ~ crpbi_hostil_rechazo + child_demo_age + child_demo_gender + (1| record_id), data = one_imp, REML =FALSE) )diag_plots_rq1 <-lapply(names(models_rq1), function(name) { df <-data.frame(res =residuals(models_rq1[[name]]),fit =fitted(models_rq1[[name]]) ) qq_rq1 <-ggplot(df, aes(sample = res)) +stat_qq(alpha =0.5) +stat_qq_line(color ="red") +labs(title = name, x ="Theoretical", y ="Sample") +theme_minimal(base_size =9) hist_rq1 <-ggplot(df, aes(x = res)) +geom_histogram(bins =20, fill ="lightblue", color ="white") +geom_vline(xintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Residuals", y ="Count") +theme_minimal(base_size =9) rvf_rq1 <-ggplot(df, aes(x = fit, y = res)) +geom_point(alpha =0.3) +geom_hline(yintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Fitted", y ="Residuals") +theme_minimal(base_size =9)list(qq = qq_rq1, rvf = rvf_rq1, hist = hist_rq1)})# q-q plots qqplots_rq1 <-wrap_plots(lapply(diag_plots_rq1, `[[`, "qq"), ncol =3)# residual variance residual_spread_rq1 <-wrap_plots(lapply(diag_plots_rq1, `[[`, "rvf"), ncol =3) hist_rq1 <-wrap_plots(lapply(diag_plots_rq1, `[[`, "hist"), ncol =3)```#### NormalityI examined the normality of model level 1 residuals visually using Q-Q plots and histograms of the residuals, and formally using the Shapiro-Wilk test. Residuals were approximately normally distributed, with no severe departures observed. Shapiro-Wilk tests indicated that residuals were normally distributed across all models (W range: 0.983–0.995), with the exception of the model predicting externalizing symptoms from communicative parenting, which showed marginal violations (W = 0.983, p = .047). Given the sensitivity of the Shapiro-Wilk test to sample size and the absence of meaningful deviations in the corresponding Q-Q plots, I did not interpret this as a substantive violation of the normality assumption.```{r}#| label: fig-qq-rq1#| fig-cap: Q-Q Plots for Examining Normality of the Level 1 Residuals for RQ1#| echo: falseqqplots_rq1``````{r}#| label: fig-histo-rq1#| fig-cap: Histograms for Examining Normality of the Level 1 Residuals for RQ1#| echo: falsehist_rq1``````{r}#| label: regression-diagnostics-rq1-l2#| echo: false#| include: false#| results: false# Level 2 EB residuals — following Cho et al. (2022)l2_plots_rq1 <-lapply(names(models_rq1), function(name) { eb <-ranef(models_rq1[[name]])$record_id[, 1] eb_std <- eb /sd(eb) df_l2 <-data.frame(eb_std = eb_std) qq_l2 <-ggplot(df_l2, aes(sample = eb_std)) +stat_qq(alpha =0.5) +stat_qq_line(color ="red") +labs(title = name, x ="Theoretical", y ="Standardized EB") +theme_minimal(base_size =9) hist_l2 <-ggplot(df_l2, aes(x = eb_std)) +geom_histogram(aes(y =after_stat(density)),bins =15, fill ="darkgrey", color ="white") +geom_vline(xintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Standardized EB (Random Intercept)", y ="Count") +theme_minimal(base_size =9)list(qq = qq_l2, hist = hist_l2)})qqplots_l2_rq1 <-wrap_plots(lapply(l2_plots_rq1, `[[`, "qq"), ncol =3)histplots_l2_rq1 <-wrap_plots(lapply(l2_plots_rq1, `[[`, "hist"), ncol =3)# Shapiro-Wilk for Level 2shapiro_l2_rq1 <-data.frame(Model =names(models_rq1), W =NA, p =NA)for (i inseq_along(models_rq1)) { eb <-ranef(models_rq1[[i]])$record_id[, 1] eb_std <- eb /sd(eb) sw <-shapiro.test(eb_std) shapiro_l2_rq1$W[i] <-round(sw$statistic, 3) shapiro_l2_rq1$p[i] <-round(sw$p.value, 3)}```I additionally examined the normality of Level 2 random effects (empirical Bayes residuals for the random intercept) using Q-Q plots, histograms, and Shapiro-Wilk tests, following @cho2022's recommendations for multilevel models.Random intercepts were approximately normally distributed across internalizing models (W range: 0.988–0.992, ps \> .05). There were mild departures in normality for externalizing models in the Shapiro-Wilk tests (W range: 0.968–0.976, ps = .015–.061), but my visual inspection of histograms revealed only slight skew without extreme deviations.Given that the data only had 100 clusters available at Level 2, the Shapiro-Wilk test has limited reliability and @cho2022 recommends visual inspection as the more informative criterion with a small number of clusters. Because of these considerations, I did not interpret these as substantive violations.```{r}#| label: fig-qq-l2-rq1#| fig-cap: "Q-Q Plots of Level 2 Empirical Bayes Residuals (Random Intercepts)"#| echo: falseqqplots_l2_rq1``````{r}#| label: fig-histo-l2-rq1#| fig-cap: "Histograms of Level 2 Empirical Bayes Residuals (Random Intercepts)"#| echo: falsehistplots_l2_rq1``````{r}#| label: shapiro-l2-rq1#| echo: falseshapiro_l2_rq1```#### HomoskedasticityI examined homoskedasticity at Level 1 only, as all models included only a random intercept at Level 2 (i.e., no Level 2 predictors. As can be seen in @fig-residuals-rq1 assumption that residuals have constant variance looked reasonably met across all 6 models, with residuals scattered fairly randomly around zero without any obvious funneling or systematic patterns.```{r}#| label: fig-residuals-rq1#| fig-cap: Examining the Distribution of Model Residuals to Assess Homoskedasticityresidual_spread_rq1```#### MulticollinearityI examined bivariate correlations among the predictor variables, computed variance inflation factors (VIF) and summed the reciprocal of the eigenvalues for all predictors included in Models 2. As seen in the heatmap below, the correlations among predictors were low to moderate (all \|r\| \< .40). VIF values were all below the conventional threshold of 5 (range: 1.01–1.16), as can be seen in the output below.```{r}#| label: multicollinearity-rq1-computations#| results: false# corrscor_vars_rq1 <- one_imp %>%select(crpbi_comunicativo, crpbi_controlador, crpbi_hostil_rechazo, child_demo_age, child_demo_gender) %>%mutate(crpbi_comunicativo =as.numeric(crpbi_comunicativo),crpbi_controlador =as.numeric(crpbi_controlador),crpbi_hostil_rechazo =as.numeric(crpbi_hostil_rechazo),child_demo_age =as.numeric(child_demo_age),child_demo_gender =as.factor(child_demo_gender) )het_cor_rq1 <-hetcor(cor_vars_rq1, use ="pairwise.complete.obs")cor_mat_rq1 <- het_cor_rq1$correlationsrownames(cor_mat_rq1) <-colnames(cor_mat_rq1) <-c("Warm/Communicative Parenting", "Controlling Parenting", "Hostile/Rejecting Parenting","Child Age", "Child Sex")# VIFscov_models_rq1 <-list("Int-Warm/Communicative + Covs"=lm(sdq_internalizing ~ crpbi_comunicativo + child_demo_age + child_demo_gender, data = one_imp),"Int-Controlling + Covs"=lm(sdq_internalizing ~ crpbi_controlador + child_demo_age + child_demo_gender, data = one_imp),"Int-Hostile + Covs"=lm(sdq_internalizing ~ crpbi_hostil_rechazo + child_demo_age + child_demo_gender, data = one_imp),"Ext-Warm/Communicative + Covs"=lm(sdq_externalizing ~ crpbi_comunicativo + child_demo_age + child_demo_gender, data = one_imp),"Ext-Controlling + Covs"=lm(sdq_externalizing ~ crpbi_controlador + child_demo_age + child_demo_gender, data = one_imp),"Ext-Hostile + Covs"=lm(sdq_externalizing ~ crpbi_hostil_rechazo + child_demo_age + child_demo_gender, data = one_imp) )vif_results <-do.call(rbind, lapply(names(cov_models_rq1), function(name) { v <-vif(cov_models_rq1[[name]])data.frame(Model = name, Predictor =names(v), VIF =round(v, 2)) }))``````{r}#| label: corr-matrix-rq1#| echo: false#| corrplot(cor_mat_rq1, method ="color", type ="upper",addCoef.col ="black", number.cex =0.7,tl.cex =0.8, tl.col ="black")``````{r}#| label: multicollinearity-rq1-results#| echo: falsevif_results``````{r}#| label: eigenvalues-rq1#| results: false### eigenvalues # getting eigenvalues of the design matrix cor_vars_rq1 <- one_imp %>%select(crpbi_comunicativo, crpbi_controlador, crpbi_hostil_rechazo, child_demo_age, child_demo_gender) %>%mutate(crpbi_comunicativo =as.numeric(crpbi_comunicativo),crpbi_controlador =as.numeric(crpbi_controlador),crpbi_hostil_rechazo =as.numeric(crpbi_hostil_rechazo),child_demo_age =as.numeric(child_demo_age),child_demo_gender =as.factor(child_demo_gender) ) het_cor_rq1 <-hetcor(cor_vars_rq1, use ="pairwise.complete.obs") cor_X_rq1 <- het_cor_rq1$correlations eigenvalues_rq1 <-eigen(cor_X_rq1)$values# summing reciprocals sum_reciprocals_rq1 <-sum(1/ eigenvalues_rq1)```Finally, I examined the sum of the reciprocals of the eigenvalues of the matrix of predictors for this second research question, and the total was `r round(sum_reciprocals_rq1, 2)`, which provided additional evidence in support of there being no collinearity issues.### Results```{r}#| label: rq1-models#| results: false###### RQ1: ICCs ####### Null models - no predictorsnull_model_int <-with(implist_std, lmer(sdq_internalizing ~1+ (1|record_id)))null_model_ext <-with(implist_std, lmer(sdq_externalizing ~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_internalizing ~ crpbi_comunicativo + (1| record_id), REML =FALSE)) model_int_comunicativo2 <-with(implist_std, lmer(sdq_internalizing ~ 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_internalizing ~ crpbi_controlador + (1| record_id), REML =FALSE)) model_int_controlador2 <-with(implist_std, lmer(sdq_internalizing ~ 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_internalizing ~ crpbi_hostil_rechazo + (1| record_id), REML =FALSE)) model_int_rechazo2 <-with(implist_std, lmer(sdq_internalizing ~ 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_externalizing ~ crpbi_comunicativo + (1| record_id), REML =FALSE)) model_ext_comunicativo2 <-with(implist_std, lmer(sdq_externalizing ~ 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_externalizing ~ crpbi_controlador + (1| record_id), REML =FALSE)) model_ext_controlador2 <-with(implist_std, lmer(sdq_externalizing ~ 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_externalizing ~ crpbi_hostil_rechazo + (1| record_id), REML =FALSE)) model_ext_rechazo2 <-with(implist_std, lmer(sdq_externalizing ~ 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)``````{r}#| label: tbl-rq2#| echo: false#| result: true#### extract fixed effects for one modelextract_fixed <-function(model, model_name, main_predictor) { est <-testEstimates(model, extra.pars =TRUE)$estimates est <-as.data.frame(est) est$term <-rownames(est) pval_col <-grep("^P", names(est), value =TRUE)[1] est <- est[!is.na(est[[pval_col]]), ] est$term[est$term == main_predictor] <-"parenting_subscale" est$CI_low <- est$Estimate -1.96* est$Std.Error est$CI_high <- est$Estimate +1.96* est$Std.Error est$b_ci <-sprintf("%.3f [%.3f, %.3f]", est$Estimate, est$CI_low, est$CI_high) est$se <-sprintf("%.3f", est$Std.Error) est$p <-ifelse(est[[pval_col]] < .001, "< .001", sprintf("%.3f", est[[pval_col]])) est$sig <- est[[pval_col]] < .05 est %>%select(term, b_ci, se, p, sig) %>%rename_with(~paste0(model_name, "_", .), -term)}#### Extract variance components for one modelextract_vc <-function(model, model_name) { extra <-as.data.frame(testEstimates(model, extra.pars =TRUE)$extra.pars) extra$term <-rownames(extra) tau00 <- extra$Estimate[grep("Intercept~~Intercept", extra$term)] sigma2 <- extra$Estimate[grep("Residual~~Residual", extra$term)] icc <- tau00 / (tau00 + sigma2)data.frame(term =c("tau00", "sigma2", "icc"),b_ci =c(sprintf("%.3f", tau00), sprintf("%.3f", sigma2), sprintf("%.3f", icc)),se =c("—", "—", "—"),p =c("—", "—", "—"),sig =c(FALSE, FALSE, FALSE) ) %>%rename_with(~paste0(model_name, "_", .), -term)}#### vc labelsvc_labels <-c("tau00"="τ00 (between-cluster)","sigma2"="σ² (within-cluster)","icc"="ICC")#### fx to build table: one per parenting subscalebuild_table_combined <-function( m_int1, m_int2, m_ext1, m_ext2, main_predictor, predictor_label, table_number) {##### Build fe_labels dynamically using the actual subscale name fe_labels <-c("(Intercept)"="Intercept","parenting_subscale"= predictor_label,"child_demo_age"="Child Age","child_demo_gender1"="Child Sex" )# Fixed effects fe_int1 <-extract_fixed(m_int1, "INT_M1", main_predictor) fe_int2 <-extract_fixed(m_int2, "INT_M2", main_predictor) fe_ext1 <-extract_fixed(m_ext1, "EXT_M1", main_predictor) fe_ext2 <-extract_fixed(m_ext2, "EXT_M2", main_predictor) fe_tbl <- fe_int1 %>%full_join(fe_int2, by ="term") %>%full_join(fe_ext1, by ="term") %>%full_join(fe_ext2, by ="term") %>%mutate(term = dplyr::recode(term, !!!fe_labels)) %>%mutate(term =factor(term, levels = fe_labels)) %>%arrange(term) %>%mutate(term =as.character(term))# Extract sig rows (BEFORE removing sig columns) sig_rows_int1 <-which(fe_tbl$INT_M1_sig ==TRUE) sig_rows_int2 <-which(fe_tbl$INT_M2_sig ==TRUE) sig_rows_ext1 <-which(fe_tbl$EXT_M1_sig ==TRUE) sig_rows_ext2 <-which(fe_tbl$EXT_M2_sig ==TRUE) parenting_row <-which(fe_tbl$term == predictor_label) stub_bold_rows <-setdiff(unique(c(sig_rows_int1, sig_rows_int2, sig_rows_ext1, sig_rows_ext2)), parenting_row ) fe_tbl <- fe_tbl %>%select(-ends_with("_sig"))#### vc_int1 <-extract_vc(m_int1, "INT_M1") vc_int2 <-extract_vc(m_int2, "INT_M2") vc_ext1 <-extract_vc(m_ext1, "EXT_M1") vc_ext2 <-extract_vc(m_ext2, "EXT_M2") vc_tbl <- vc_int1 %>%full_join(vc_int2, by ="term") %>%full_join(vc_ext1, by ="term") %>%full_join(vc_ext2, by ="term") %>%mutate(term = dplyr::recode(term, !!!vc_labels)) %>%mutate(term =factor(term, levels = vc_labels)) %>%arrange(term) %>%mutate(term =as.character(term)) %>%select(-ends_with("_sig"))#### combining and replacing missing vals n_fe <-nrow(fe_tbl) n_vc <-nrow(vc_tbl) full_tbl <-bind_rows( fe_tbl %>%mutate(row_group ="Fixed Effects"), vc_tbl %>%mutate(row_group ="Variance Components") ) %>%mutate(across(-c(term, row_group), ~ifelse(is.na(.), "—", .)))## gt tbl full_tbl %>%select(-row_group) %>%gt(rowname_col ="term") %>%# Row groupstab_row_group(label ="Variance Components", rows = (n_fe +1):(n_fe + n_vc)) %>%tab_row_group(label ="Fixed Effects", rows =1:n_fe) %>%row_group_order(groups =c("Fixed Effects", "Variance Components")) %>%# Top-level spanners — outcometab_spanner(label ="Internalizing Symptoms", id ="spanner_int", columns =starts_with("INT")) %>%tab_spanner(label ="Externalizing Symptoms", id ="spanner_ext", columns =starts_with("EXT")) %>%# Nested model spannerstab_spanner(label ="Model 1", id ="int_m1", columns =c(INT_M1_b_ci, INT_M1_se, INT_M1_p)) %>%tab_spanner(label ="Model 2", id ="int_m2", columns =c(INT_M2_b_ci, INT_M2_se, INT_M2_p)) %>%tab_spanner(label ="Model 1", id ="ext_m1", columns =c(EXT_M1_b_ci, EXT_M1_se, EXT_M1_p)) %>%tab_spanner(label ="Model 2", id ="ext_m2", columns =c(EXT_M2_b_ci, EXT_M2_se, EXT_M2_p)) %>%cols_label(INT_M1_b_ci ="β [95% CI]", INT_M1_se ="SE", INT_M1_p ="p",INT_M2_b_ci ="β [95% CI]", INT_M2_se ="SE", INT_M2_p ="p",EXT_M1_b_ci ="β [95% CI]", EXT_M1_se ="SE", EXT_M1_p ="p",EXT_M2_b_ci ="β [95% CI]", EXT_M2_se ="SE", EXT_M2_p ="p" ) %>%tab_stubhead(label ="Predictor") %>%tab_source_note(source_note ="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)." ) %>%tab_options(table.border.top.style ="solid",table.border.top.width =px(2),table.border.bottom.style ="solid",table.border.bottom.width =px(2),heading.border.bottom.style ="solid",heading.border.bottom.width =px(1),column_labels.border.top.style ="none",column_labels.border.bottom.style ="solid",column_labels.border.bottom.width =px(1),table_body.border.top.style ="none",table_body.border.bottom.style ="none",stub.border.style ="none",row_group.border.top.style ="solid",row_group.border.top.width =px(1),row_group.border.bottom.style ="none",row.striping.include_table_body =FALSE,table.font.size =px(12),table.font.names ="Times New Roman",heading.title.font.size =px(12),heading.subtitle.font.size =px(12),heading.align ="left",data_row.padding =px(4) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_title("title") ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_title("subtitle") ) %>%tab_style(style =cell_text(align ="left"),locations =cells_stub() ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_row_groups() ) %>%#### bolding significancetab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(INT_M1_b_ci, INT_M1_p), rows = sig_rows_int1) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(INT_M2_b_ci, INT_M2_p), rows = sig_rows_int2) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(EXT_M1_b_ci, EXT_M1_p), rows = sig_rows_ext1) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(EXT_M2_b_ci, EXT_M2_p), rows = sig_rows_ext2) )}#### tablestbl_comunicativo <-build_table_combined(m_int1 = model_int_comunicativo1, m_int2 = model_int_comunicativo2,m_ext1 = model_ext_comunicativo1, m_ext2 = model_ext_comunicativo2,main_predictor ="crpbi_comunicativo",predictor_label ="Communicative Parenting",table_number =9)tbl_controlador <-build_table_combined(m_int1 = model_int_controlador1, m_int2 = model_int_controlador2,m_ext1 = model_ext_controlador1, m_ext2 = model_ext_controlador2,main_predictor ="crpbi_controlador",predictor_label ="Controlling Parenting",table_number =10)tbl_rechazo <-build_table_combined(m_int1 = model_int_rechazo1, m_int2 = model_int_rechazo2,m_ext1 = model_ext_rechazo1, m_ext2 = model_ext_rechazo2,main_predictor ="crpbi_hostil_rechazo",predictor_label ="Hostile/Rejecting Parenting",table_number =11)```#### Warm/Communicative ParentingModel 1 examined the bivariate association between communicative parenting and each outcome, controlling for family clustering. Warm/communicative parenting was not a significant predictor of internalizing symptoms in Model 1 (β = -0.079, SE = 0.078, p = .314), nor did it reach significance in Model 2 after adding child age and child sex (β = -0.092, SE = 0.075, p = .219).For externalizing symptoms, communicative parenting was a significant negative predictor in Model 1 (β = -0.169, SE = 0.078, p = .031), such that youth reporting higher levels of warm/communicative parenting showed lower externalizing symptoms. This effect remained significant and stable in Model 2 (β = -0.171, SE = 0.079, p = .030).In terms of child covariates, both youth age and sex emerged as significant negative predictors of internalizing symptoms (β = -0.229, SE = 0.074, p = .002; β = 0.382, SE = 0.145, p = .008). This indicated that older and male youth reported lower internalizing symptoms.For externalizing symptoms neither youth covariate was significant. However, child sex showed a near-significant trend (β = -0.268, SE = 0.153, p = .079), suggesting the possibility that female youth may report lower externalizing symptoms in this sample.These results provide partial support for Hypothesis 1 — warm/communicative parenting was associated with lower externalizing symptoms as predicted, but the association with internalizing symptoms was not significant.```{r}#| label: tbl-com-int-ext#| tbl-cap: "Pooled Multilevel Model Estimates for Warm/Communicative Parenting as a Predictor of Internalizing and Externalizing Behavior"tbl_comunicativo```#### Controlling ParentingControlling parenting was a significant positive predictor of both internalizing (β = 0.325, SE = 0.075, p \< .001) and externalizing symptoms (β = 0.182, SE = 0.078, p = .020), such that youth reporting higher levels of controlling parenting showed higher levels of both internalizing and externalizing symptoms.In Model 2, after adding child age and child sex, the effect of controlling parenting strengthened slightly for internalizing (β = 0.384, SE = 0.076, p \< .001) and remained stable for externalizing (β = 0.191, SE = 0.083, p = .021).Child sex emerged as a significant predictor of internalizing symptoms in Model 2 (β = 0.601, SE = 0.139, p \< .001), indicating that female youth reported higher internalizing symptoms when accounting for controlling parenting and age. Child age showed a near-significant negative trend for internalizing (β = -0.129, SE = 0.071, p = .069). These results provide support for Hypothesis 3 — controlling parenting was associated with both internalizing and externalizing symptoms. In this sample, the direction of the association was positive across both outcomes.```{r}#| label: tbl-con-int-ext#| tbl-cap: "Pooled Multilevel Model Estimates for Controlling Parenting as a Predictor of Internalizing and Externalizing Behavior"tbl_controlador```#### Hostile/Rejecting Parenting```{r}#| label: tbl-hos-int-ext#| tbl-cap: "Pooled Multilevel Model Estimates for Hostile/Rejecting Parenting as a Predictor of Internalizing and Externalizing Behavior"tbl_rechazo```Hostile/rejecting parenting was a significant positive predictor of both internalizing (β = 0.248, SE = 0.076, p = .001) and externalizing symptoms (β = 0.393, SE = 0.073, p \< .001), such that youth reporting higher levels of hostile and rejecting parenting showed higher levels of both internalizing and externalizing symptoms.In Model 2, after adding child age and child sex, the effect of hostile/rejecting parenting remained stable and significant for both internalizing (β = 0.256, SE = 0.072, p \< .001) and externalizing (β = 0.395, SE = 0.073, p \< .001).Child age emerged as a significant negative predictor of internalizing symptoms in Model 2 (β = -0.192, SE = 0.072, p = .007), with older youth reporting lower internalizing symptoms. Child sex was also significant for internalizing (β = 0.459, SE = 0.140, p = .001), indicating that female youth reported higher internalizing symptoms when accounting for hostile/rejecting parenting and age. Neither child age nor child sex significantly predicted externalizing symptoms in Model 2.These results provide support for Hypothesis 2 — hostile and rejecting parenting was associated with higher internalizing and externalizing symptoms in this sample, as predicted.#### Overall FindingsAcross all three parenting dimensions, the pattern of results provided at least partial support for the hypotheses. Hypothesis 1, which predicted that higher levels of parental warmth would be associated with lower internalizing and externalizing symptoms, was partially supported — warm/communicative parenting significantly predicted lower externalizing symptoms but was not a significant predictor of internalizing symptoms.Hypothesis 2, which predicted that hostile and rejecting parenting would be associated with higher internalizing and externalizing symptoms, was fully supported — hostile/rejecting parenting was a significant positive predictor of both outcomes across both models. Hypothesis 3, which proposed that controlling parenting would be associated with youth symptoms with the direction to be determined empirically, was also supported — controlling parenting was a significant predictor of both internalizing and externalizing symptoms across both models. Importantly, this was a positive association so that higher levels of controlling parenting predicted higher levels of both internalizing and externalizing symptoms. These effects were of higher magnitude once child age and sex were included in the models.Notably, hostile/rejecting parenting had a stronger effect on youth externalizing in comparison to internalizing symptoms, while the effect of controlling parenting was of higher magnitude for youth internalizing symptoms.Child demographic characteristics emerged as consistent covariates for internalizing but not externalizing symptoms. Specifically, female and younger youth were more likely to report higher internalizing symptoms, once parenting dimensions were accounted for. The ICC values across models ranged from 0.040 to 0.280, indicating meaningful family-level clustering in youth outcomes, particularly for internalizing symptoms.Model comparisons using the D3 pooled likelihood ratio test (@meng1992) indicated that adding child age and sex in Model 2 significantly improved model fit for internalizing symptoms across all three parenting dimensions: warm/communicative parenting (F(2) = 8.23, p \< .001), controlling parenting (F(2) = 11.15, p \< .001), and hostile/rejecting parenting (F(2) = 8.88, p \< .001). In contrast, the addition of child age and sex did not significantly improve model fit for externalizing symptoms across any parenting dimension: warm/communicative parenting (F(2) = 1.82, p = .162), controlling parenting (F(2) = 1.68, p = .187), and hostile/rejecting parenting (F(2) = 1.90, p = .150) — suggesting that child demographic characteristics accounted for unique variance in internalizing but not externalizing symptoms beyond parenting.## RQ2: Identifying Risk factors and Promotive Factors for Mothers’ Parental Resilience in Informal SettlementsI drew on @gavidiapayne2015 to define parental resilience as “the capacity of parents to deliver competent, quality parenting to children despite adverse personal, family, and social circumstances.” (p. 113). Specifically, parental resilience regarding positive youth mental health outcomes. Drawing on RQ1 results, this meant high warm/communicative parenting, and low controlling and hostile/rejecting parenting.I examined individual, relational, and community-level factors that either hindered (risk factors) or supported (promotive factors) mothers’ parental resilience in this context.**Hypothesis 4:** Individual-level factors will predict parental resilience (good parenting quality) after controlling for youth characteristics and family demographics. Specifically, lower maternal mental health symptoms, higher educational attainment, age, will be associated with higher parental resilience (higher warmth, lower hostility), with the direction for controlling parenting to be determined empirically.**Hypothesis 5:** Relational-level factors will predict parental resilience above and beyond individual-level factors and covariates. Specifically, higher social connections will be associated with higher parental resilience (higher warmth, lower hostility), with the direction for controlling parenting to be determined empirically. This tests whether relational resources explain unique variance in parenting quality beyond individual characteristics and adversity exposure.**Hypothesis 6:** Community-level factors will predict parental resilience above and beyond individual and relational factors. Specifically, higher perceived community mutual aid and/or participation will be associated with higher parental resilience (higher warmth, lower hostility), with the direction for controlling parenting to be determined empirically. This tests whether community-level resources contribute to parenting quality beyond individual and relational determinants.**Hypothesis 7:** According to resilience frameworks that consider both risks and protective factors, every ecological level is expected to make a unique contribution to parenting quality. More specifically, relational factors (Model 3) and community factors (Model 4) should add significant predictive value for parental resilience beyond what individual factors alone (Model 2) can explain. This will be shown through likelihood ratio tests, and greater reductions in Level 2 variance.Following these hypotheses, four sequential models tested unique variance explained by each ecological level. The goal with this approach was to examine whether more distal factors (e.g., 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.All Level 2 predictors were grand mean centered and the models were estimated employing Maximum Likelihood Estimation (MLE) using the lme4 package in R (@bates2015) 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 (@finch2019). I used these final estimates to compare with the standard OLS estimates in the final step of analyses. Formal model equations follow.::: {style="text-align: center"}```{=latex}\usepackage{booktabs}\begin{table}[htbp]\caption*{\textit{Proposed Models for Research Question 2}}\centering\renewcommand{\arraystretch}{2.0}\begin{tabular}{ccp{11cm}r}\toprule\textbf{Model} & \textbf{Level} & \multicolumn{1}{c}{\textbf{Equation}} & \\\midrule0 & 1 & $\mathit{PQual}_{ij} = \beta_{0j} + r_{ij}$ & (Eq. 1) \\ & 2 & $\beta_{0j} = \gamma_{00} + u_{0j}$ & (Eq. 2) \\[0.5em]1 & 1 & $\mathit{PQuality}_{ij} = \beta_{0j} + \beta_1(\mathit{TotProbs})_{ij} + \beta_2(\mathit{ChildAge})_{ij} + \beta_3(\mathit{ChildFemale})_{ij} + r_{ij}$ & (Eq. 3) \\ & 2 & $\beta_{0j} = \gamma_{00} + \gamma_{01}(\mathit{MotherEdu})_j + u_{0j}$ & (Eq. 4) \\[0.5em]2 & 1 & {[Same as Model 1]} & \\ & 2 & $\beta_{0j} = \gamma_{00} + \gamma_{01}(\mathit{MotherEdu})_j + \gamma_{02}(\mathit{MotherChildAdv})_j + \gamma_{03}(\mathit{InfSettlAdv})_j + \gamma_{04}(\mathit{MotherMH})_j + u_{0j}$ & (Eq. 5) \\[0.5em]3 & 1 & {[Same as Model 1]} & \\ & 2 & $\beta_{0j} = \gamma_{00} + \gamma_{01}(\mathit{MotherEdu})_j + \gamma_{02}(\mathit{MotherChildAdv})_j + \gamma_{03}(\mathit{InfSettlAdv})_j + \gamma_{04}(\mathit{MotherMH})_j + \gamma_{05}(\mathit{MotherSocCon})_j + u_{0j}$ & (Eq. 6) \\[0.5em]4 & 1 & {[Same as Model 1]} & \\ & 2 & $\beta_{0j} = \gamma_{00} + \gamma_{01}(\mathit{MotherEdu})_j + \gamma_{02}(\mathit{MotherChildAdv})_j + \gamma_{03}(\mathit{InfSettlAdv})_j + \gamma_{04}(\mathit{MotherMH})_j + \gamma_{05}(\mathit{MotherSocCon})_j + \gamma_{06}(\mathit{CommOrg})_j + \gamma_{07}(\mathit{CommMutualAid})_j + u_{0j}$ & (Eq. 7) \\\bottomrule\end{tabular}\vspace{1em}\begin{center}\textbf{Where:}\\[0.5em]\begin{tabular}{rl}$\mathit{PQuality}_{ij}$ & $=$ Parenting dimension score for youth $i$ in family $j$ \\$\beta_1$ & $=$ Fixed effect of youth total problems \\$\beta_2$ & $=$ Fixed effect of youth age \\$\beta_3$ & $=$ Fixed effect of youth sex (female) \\$\mathit{TotProbs}_{ij}$ & $=$ Total problems score for youth $i$ in family $j$ \\$\mathit{ChildAge}_{ij}$ & $=$ Age of youth $i$ in family $j$ \\$\mathit{ChildFemale}_{ij}$ & $=$ 1 if youth $i$ in family $j$ is female \\$r_{ij}$ & $=$ Level 1 residual (individual youth deviation from family mean) \\[0.5em]$\gamma_{00}$ & $=$ Fixed intercept (grand mean of parenting dimension across all youth) \\$\gamma_{01}$ & $=$ Fixed effect of maternal years of education \\$\gamma_{02}$ & $=$ Fixed effect of maternal childhood adversity \\$\gamma_{03}$ & $=$ Fixed effect of informal settlement adversity \\\end{tabular}\end{center}\end{table}```:::### Regression AssumptionsI conducted a series of analyses to identify any potential violations of the regression assumptions of linearity, normality, homoskedasticity, and multicollinearity.#### Linear RelationshipI examined Level 1 linearity by plotting marginal standardized residuals against each continuous Level 1 predictor, following @cho2022) (see @fig-l1-linearity-sdq-rq2 and @fig-l1-linearity-age-rq2). Residuals were generally scattered around zero across the full range of both SDQ total problems and child age for all three parenting outcomes, with no systematic nonlinear patterns observed, supporting the assumption of linearity for Level 1 predictors. A small number of observations showed marginal standardized residuals beyond ±3, most notably in the Communicative Parenting model where a few observations approached ±4. Given the bounded 1–3 scale of the parenting outcomes, I attributed these extreme residuals to the discrete distributional constraints of the outcome rather than true model violations, and noted them as a limitation of applying linear models to bounded Likert-type outcomes. I assessed Level 2 linearity via plots of standardized EB residuals against each Level 2 predictor (see @fig-l2-hetero-comm, @fig-l2-hetero-ctrl, and @fig-l2-hetero-host), which showed no systematic patterns inconsistent with linearity.```{r}#| label: rq2-l1-linearity#| echo: false#| results: falsemodels_rq2 <-list("Warm/Communicative M4"=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), data = one_imp, REML =TRUE),"Controlling M4"=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), data = one_imp, REML =TRUE),"Hostile/Rejecting M4"=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), data = one_imp, REML =TRUE))l1_linearity_rq2 <-lapply(names(models_rq2), function(name) {# Marginal residuals = conditional residuals + random effects cond_res <-residuals(models_rq2[[name]], type ="pearson") ran_ef <-ranef(models_rq2[[name]])$record_id[, 1] re_vals <- ran_ef[match(one_imp$record_id,as.integer(rownames(ranef(models_rq2[[name]])$record_id)))] mar_res <- cond_res + re_vals mar_res_std <- mar_res /sd(mar_res, na.rm =TRUE) df <-data.frame(mar_res_std = mar_res_std,sdq_total = one_imp$sdq_total,child_age = one_imp$child_demo_age ) p_sdq <-ggplot(df, aes(x = sdq_total, y = mar_res_std)) +geom_point(alpha =0.3) +geom_hline(yintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="SDQ Total", y ="Marginal Std Residuals") +theme_minimal(base_size =9) p_age <-ggplot(df, aes(x = child_age, y = mar_res_std)) +geom_point(alpha =0.3) +geom_hline(yintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Child Age", y ="Marginal Std Residuals") +theme_minimal(base_size =9)list(sdq = p_sdq, age = p_age)})linearity_sdq_rq2 <-wrap_plots(lapply(l1_linearity_rq2, `[[`, "sdq"), ncol =3)linearity_age_rq2 <-wrap_plots(lapply(l1_linearity_rq2, `[[`, "age"), ncol =3)``````{r}#| label: fig-l1-linearity-sdq-rq2#| fig-cap: "Level 1 Linearity: Marginal Standardized Residuals vs SDQ Total Problems (RQ2)"#| fig-cap-location: top#| echo: falselinearity_sdq_rq2``````{r}#| label: fig-l1-linearity-age-rq2#| fig-cap: "Level 1 Linearity: Marginal Standardized Residuals vs Child Age (RQ2)"#| fig-cap-location: top#| echo: falselinearity_age_rq2``````{r}#| label: regression-diagnostics-rq2#| echo: false#| include: false#| results: false### fitting all final models with just this one datasetmodels_rq2 <-list("Warm/Communicative M4"=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), data = one_imp, REML =TRUE),"Controlling M4"=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), data = one_imp, REML =TRUE),"Hostile/Rejecting M4"=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), data = one_imp, REML =TRUE))diag_plots_rq2 <-lapply(names(models_rq2), function(name) { df <-data.frame(res =residuals(models_rq2[[name]], type ="pearson"),fit =fitted(models_rq2[[name]]) ) qq_rq2 <-ggplot(df, aes(sample = res)) +stat_qq(alpha =0.5) +stat_qq_line(color ="red") +labs(title = name, x ="Theoretical", y ="Sample") +theme_minimal(base_size =9) hist_rq2 <-ggplot(df, aes(x = res)) +geom_histogram(aes(y =after_stat(density)),bins =20, fill ="lightblue", color ="white") +geom_vline(xintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Residuals", y ="Density") +theme_minimal(base_size =9) rvf_rq2 <-ggplot(df, aes(x = fit, y = res)) +geom_point(alpha =0.3) +geom_hline(yintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Fitted", y ="Residuals") +theme_minimal(base_size =9)list(qq = qq_rq2, rvf = rvf_rq2, hist = hist_rq2)})qqplots_rq2 <-wrap_plots(lapply(diag_plots_rq2, `[[`, "qq"), ncol =3)residual_spread_rq2 <-wrap_plots(lapply(diag_plots_rq2, `[[`, "rvf"), ncol =3)hist_rq2 <-wrap_plots(lapply(diag_plots_rq2, `[[`, "hist"), ncol =3)# Shapiro-Wilk Level 1shapiro_rq2 <-data.frame(Model =names(models_rq2), W =NA, p =NA)for (i inseq_along(models_rq2)) { res_std <-residuals(models_rq2[[i]], type ="pearson") sw <-shapiro.test(res_std) shapiro_rq2$W[i] <-round(sw$statistic, 3) shapiro_rq2$p[i] <-round(sw$p.value, 3)}```#### NormalityI examined the normality of model Level 1 residuals visually using Q-Q plots and histograms of the residuals, and formally using the Shapiro-Wilk test. Shapiro-Wilk tests indicated significant departures from normality for the Warm/Communicative (W = 0.929, p \< .001) and Hostile/Rejecting models (W = 0.974, p = .004), while the Controlling model showed no significant violation (W = 0.986, p = .103). Visual inspection of the histograms revealed left skew across all three models, most pronounced for Warm/Communicative parenting, which is consistent with the ceiling effect observed in the descriptive statistics for this variable. Given that these departures are attributable to the bounded 1–3 scale of the parenting outcomes rather than model misspecification, and that multilevel models are generally robust to moderate violations of normality with respect to fixed effect estimates (as discussed by @cho2022), I did not interpret these as substantive violations.```{r}#| label: fig-qq-rq2#| fig-cap: "Q-Q Plots for Examining Normality of the Level 1 Residuals for RQ2"#| echo: falseqqplots_rq2``````{r}#| label: fig-histo-rq2#| fig-cap: "Histograms for Examining Normality of the Level 1 Residuals for RQ2"#| echo: falsehist_rq2``````{r}#| label: shapiro-wilk-l1#| echo: falseshapiro_rq2``````{r}#| label: regression-diagnostics-rq2-l2#| echo: false#| include: false#| results: false# Level 2 EB residuals — following Cho et al. (2022)l2_plots_rq2 <-lapply(names(models_rq2), function(name) { eb <-ranef(models_rq2[[name]])$record_id[, 1] eb_std <- eb /sd(eb) df_l2 <-data.frame(eb_std = eb_std) qq_l2 <-ggplot(df_l2, aes(sample = eb_std)) +stat_qq(alpha =0.5) +stat_qq_line(color ="red") +labs(title = name, x ="Theoretical", y ="Standardized EB") +theme_minimal(base_size =9) hist_l2 <-ggplot(df_l2, aes(x = eb_std)) +geom_histogram(bins =15, fill ="darkgrey", color ="white") +geom_vline(xintercept =0, color ="red", linetype ="dashed") +labs(title = name, x ="Standardized EB (Random Intercept)", y ="Count") +theme_minimal(base_size =9)list(qq = qq_l2, hist = hist_l2)})qqplots_l2_rq2 <-wrap_plots(lapply(l2_plots_rq2, `[[`, "qq"), ncol =3)histplots_l2_rq2 <-wrap_plots(lapply(l2_plots_rq2, `[[`, "hist"), ncol =3)# Shapiro-Wilk Level 2shapiro_l2_rq2 <-data.frame(Model =names(models_rq2), W =NA, p =NA)for (i inseq_along(models_rq2)) { eb <-ranef(models_rq2[[i]])$record_id[, 1] eb_std <- eb /sd(eb) sw <-shapiro.test(eb_std) shapiro_l2_rq2$W[i] <-round(sw$statistic, 3) shapiro_l2_rq2$p[i] <-round(sw$p.value, 3)}```I additionally examined the normality of Level 2 random effects (empirical Bayes residuals for the random intercept) using Q-Q plots, histograms, and Shapiro-Wilk tests, following @cho2022's recommendations for multilevel models. Shapiro-Wilk tests indicated that residuals were normally distributed across all models (W range: .977–.993), with no models showing significant violations (all ps \> .05).```{r}#| label: fig-qq-l2-rq2#| fig-cap: "Q-Q Plots of Level 2 Empirical Bayes Residuals (Random Intercepts) for RQ2"#| echo: falseqqplots_l2_rq2``````{r}#| label: fig-histo-l2-rq2#| fig-cap: "Histograms of Level 2 Empirical Bayes Residuals (Random Intercepts) for RQ2"#| echo: falsehistplots_l2_rq2``````{r}#| label: shapiro-l2-rq2#| echo: falseshapiro_l2_rq2```#### HomoskedasticityI examined homoskedasticity at both Level 1 and Level 2, as models included Level 2 predictors. As can be seen in @fig-residuals-rq2, the assumption that residuals have constant variance looked reasonably met across all three models, with residuals scattered fairly randomly around zero without any obvious funneling or systematic patterns. At Level 2, I examined whether the variance of the standardized EB residuals was constant across the range of each Level 2 predictor (see @fig-l2-hetero-comm, @fig-l2-hetero-ctrl, and @fig-l2-hetero-host). The Hostile/Rejecting model showed no substantive concerns. The Communicative Parenting model showed mild curvature for some predictors, most evident at the extremes of the distributions where data are sparse. The Controlling Parenting model showed the most pronounced patterns, particularly for perceived mutual aid and social connections. These patterns are noted as a limitation in the interpretation of findings for controlling parenting specifically.```{r}#| label: fig-residuals-rq2#| fig-cap: "Examining the Distribution of Model Residuals to Assess Level 1 Homoskedasticity for RQ2"#| echo: falseresidual_spread_rq2``````{r}#| label: rq2-l2-hetero-computations#| echo: false#| include: false#| results: falsel2_predictors <-c("adult_education", "mother_aces", "informal_settlement_index","mother_dsmtotal", "mother_sc", "commparticipation", "mutual_aid")l2_pred_labels <-c("Maternal Education", "Maternal ACEs", "Informal Settlement Index","DSM Cross-Cutting", "Social Connections","Community Participation", "Perceived Mutual Aid")l2_hetero_plots_rq2 <-lapply(names(models_rq2), function(name) { eb <-ranef(models_rq2[[name]])$record_id[, 1] eb_std <- eb /sd(eb) mother_data <- one_imp %>%group_by(record_id) %>%slice(1) %>%ungroup() eb_df <-data.frame(record_id =as.integer(rownames(ranef(models_rq2[[name]])$record_id)),eb_std = eb_std ) plot_df <-left_join(eb_df, mother_data, by ="record_id") plots <-lapply(seq_along(l2_predictors), function(j) { pred <- l2_predictors[j] label <- l2_pred_labels[j]ggplot(plot_df, aes(x = .data[[pred]], y = eb_std)) +geom_point(alpha =0.4) +geom_hline(yintercept =0, color ="red", linetype ="dashed") +labs(title = label, x = label, y ="Std EB") +theme_minimal(base_size =8) })wrap_plots(plots, ncol =4)})l2_hetero_comm <- l2_hetero_plots_rq2[[1]]l2_hetero_ctrl <- l2_hetero_plots_rq2[[2]]l2_hetero_host <- l2_hetero_plots_rq2[[3]]``````{r}#| label: fig-l2-hetero-comm#| fig-cap: "Level 2 Heteroskedasticity: Standardized EB Residuals vs Level 2 Predictors for Communicative Parenting"#| echo: falsel2_hetero_comm``````{r}#| label: fig-l2-hetero-ctrl#| fig-cap: "Level 2 Heteroskedasticity: Standardized EB Residuals vs Level 2 Predictors for Controlling Parenting"#| echo: falsel2_hetero_ctrl``````{r}#| label: fig-l2-hetero-host#| fig-cap: "Level 2 Heteroskedasticity: Standardized EB Residuals vs Level 2 Predictors for Hostile/Rejecting Parenting"#| echo: falsel2_hetero_host```#### MulticollinearityI examined bivariate correlations among the predictor variables, computed variance inflation factors (VIF), and summed the reciprocal of the eigenvalues for all predictors included in Model 4. As seen in the heatmap below, the correlations among predictors were low to moderate (all \|r\| \< .50). VIF values were all below the conventional threshold of 5 (range: 1.05–1.29), as can be seen in the output below.```{r}#| label: multicollinearity-rq2-computations#| results: false#| cor_vars_rq2 <- one_imp %>%group_by(record_id) %>%slice(1) %>%ungroup() %>%select(sdq_total, child_demo_age, child_demo_gender, adult_education, mother_aces, informal_settlement_index, mother_dsmtotal, mother_sc, commparticipation, mutual_aid) %>%mutate(across(everything(), ~as.numeric(as.character(.))))cor_mat_rq2 <-cor(cor_vars_rq2, use ="pairwise.complete.obs")rownames(cor_mat_rq2) <-colnames(cor_mat_rq2) <-c("SDQ Total", "Child Age", "Child Sex","Maternal Education", "Maternal ACEs", "Informal Settlement","DSM Cross-Cutting", "Social Connections","Community Participation", "Mutual Aid")cor_X_rq2 <- cor_mat_rq2eigenvalues_rq2 <-eigen(cor_X_rq2)$valuessum_reciprocals_rq2 <-sum(1/ eigenvalues_rq2)cov_models_rq2 <-list("Warm/Communicative M4"=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 = one_imp),"Controlling M4"=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 = one_imp),"Hostile/Rejecting M4"=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 = one_imp))vif_results_rq2 <-do.call(rbind, lapply(names(cov_models_rq2), function(name) { v <-vif(cov_models_rq2[[name]])data.frame(Model = name, Predictor =names(v), VIF =round(v, 2))}))``````{r}#| label: corr-matrix-rq2#| echo: falsecorrplot(cor_mat_rq2, method ="color", type ="upper",addCoef.col ="black", number.cex =0.6,tl.cex =0.7, tl.col ="black")``````{r}#| label: multicollinearity-rq2-results#| echo: falsevif_results_rq2``````{r}#| label: eigenvalues-rq2#| results: falsesum_reciprocals_rq2```Finally, I examined the sum of the reciprocals of the eigenvalues of the matrix of predictors for this research question, and the total was `r round(sum_reciprocals_rq2, 2)`, which provided additional evidence in support of there being no collinearity issues.```{r}#| echo: true#| results: false###### 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 estimatestestEstimates(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 estimatestestEstimates(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 estimatestestEstimates(hos_m1, extra.pars =TRUE)testEstimates(hos_m2, extra.pars =TRUE)testEstimates(hos_m3, extra.pars =TRUE)testEstimates(hos_m4_ml, extra.pars =TRUE)``````{r}#| echo: true#| results: false 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))``````{r}#| label: creating-tables#| echo: false#| include: false###function to extract pooled estimatesextract_estimates_2 <-function(model, model_name) { est <-testEstimates(model, extra.pars =TRUE)$estimates est <-as.data.frame(est) est$term <-rownames(est) pval_col <-grep("^P", names(est), value =TRUE)[1] est <- est[!is.na(est[[pval_col]]), ] est$CI_low <- est$Estimate -1.96* est$Std.Error est$CI_high <- est$Estimate +1.96* est$Std.Error est$b_ci <-sprintf("%.3f [%.3f, %.3f]", est$Estimate, est$CI_low, est$CI_high) est$se <-sprintf("%.3f", est$Std.Error) est$p <-ifelse(est[[pval_col]] < .001, "< .001", sprintf("%.3f", est[[pval_col]])) est$sig <- est[[pval_col]] < .05 est %>%select(term, b_ci, se, p, sig) %>%rename_with(~paste0(model_name, "_", .), -term)}#### Extract variance componentsextract_vc <-function(model, model_name) { extra <-as.data.frame(testEstimates(model, extra.pars =TRUE)$extra.pars) extra$term <-rownames(extra) tau00 <- extra$Estimate[grep("Intercept~~Intercept", extra$term)] sigma2 <- extra$Estimate[grep("Residual~~Residual", extra$term)] icc <- tau00 / (tau00 + sigma2)data.frame(term =c("tau00", "sigma2", "icc"),b_ci =c(sprintf("%.3f", tau00), sprintf("%.3f", sigma2), sprintf("%.3f", icc)),se =c("—", "—", "—"),p =c("—", "—", "—"),sig =c(FALSE, FALSE, FALSE) ) %>%rename_with(~paste0(model_name, "_", .), -term)}# cleann predictor labelspredictor_labels_rq2 <-c("(Intercept)"="Intercept","sdq_total"="SDQ Total","child_demo_age"="Child Age","child_demo_gender1"="Child Sex","adult_education"="Maternal Education","mother_aces"="Maternal ACEs","informal_settlement_index"="Informal Settlement Index","mother_dsmtotal"="DSM Cross-Cutting Symptom","mother_sc"="Mother Social Connections","commparticipation"="Community Participation","mutual_aid"="Perceived Mutual Aid")vc_labels <-c("tau00"="τ00 (between-cluster)","sigma2"="σ² (within-cluster)","icc"="ICC")#### tablesbuild_table_rq2 <-function(m1, m2, m3, m4_ml, outcome_label, table_number) { e1 <-extract_estimates_2(m1, "M1") e2 <-extract_estimates_2(m2, "M2") e3 <-extract_estimates_2(m3, "M3") e4 <-extract_estimates_2(m4_ml, "M4") fe_tbl <- e1 %>%full_join(e2, by ="term") %>%full_join(e3, by ="term") %>%full_join(e4, by ="term") %>%mutate(term = dplyr::recode(term, !!!predictor_labels_rq2)) %>%mutate(term =factor(term, levels = predictor_labels_rq2)) %>%arrange(term) %>%mutate(term =as.character(term))### significance sig_rows_m1 <-which(fe_tbl$M1_sig ==TRUE) sig_rows_m2 <-which(fe_tbl$M2_sig ==TRUE) sig_rows_m3 <-which(fe_tbl$M3_sig ==TRUE) sig_rows_m4 <-which(fe_tbl$M4_sig ==TRUE) fe_tbl <- fe_tbl %>%select(-ends_with("_sig")) %>%mutate(across(everything(), ~ifelse(is.na(.), "—", .)))#### variance vc1 <-extract_vc(m1, "M1") vc2 <-extract_vc(m2, "M2") vc3 <-extract_vc(m3, "M3") vc4 <-extract_vc(m4_ml, "M4") vc_tbl <- vc1 %>%full_join(vc2, by ="term") %>%full_join(vc3, by ="term") %>%full_join(vc4, by ="term") %>%mutate(term = dplyr::recode(term, !!!vc_labels)) %>%mutate(term =factor(term, levels = vc_labels)) %>%arrange(term) %>%mutate(term =as.character(term)) %>%select(-ends_with("_sig")) %>%mutate(across(everything(), ~ifelse(is.na(.), "—", .)))#### combining n_fe <-nrow(fe_tbl) n_vc <-nrow(vc_tbl) full_tbl <-bind_rows( fe_tbl %>%mutate(row_group ="Fixed Effects"), vc_tbl %>%mutate(row_group ="Variance Components") )##### bilding tbls full_tbl %>%select(-row_group) %>%gt(rowname_col ="term") %>%# Row groupstab_row_group(label ="Variance Components", rows = (n_fe +1):(n_fe + n_vc)) %>%tab_row_group(label ="Fixed Effects", rows =1:n_fe) %>%row_group_order(groups =c("Fixed Effects", "Variance Components")) %>%tab_spanner(label ="Model 1", columns =starts_with("M1")) %>%tab_spanner(label ="Model 2", columns =starts_with("M2")) %>%tab_spanner(label ="Model 3", columns =starts_with("M3")) %>%tab_spanner(label ="Model 4", columns =starts_with("M4")) %>%cols_label(M1_b_ci ="β [95% CI]", M1_se ="SE", M1_p ="p",M2_b_ci ="β [95% CI]", M2_se ="SE", M2_p ="p",M3_b_ci ="β [95% CI]", M3_se ="SE", M3_p ="p",M4_b_ci ="β [95% CI]", M4_se ="SE", M4_p ="p" ) %>%tab_stubhead(label ="Predictor") %>%tab_source_note(source_note ="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)." ) %>%tab_options(table.border.top.style ="solid",table.border.top.width =px(2),table.border.bottom.style ="solid",table.border.bottom.width =px(2),heading.border.bottom.style ="solid",heading.border.bottom.width =px(1),column_labels.border.top.style ="none",column_labels.border.bottom.style ="solid",column_labels.border.bottom.width =px(1),table_body.border.top.style ="none",table_body.border.bottom.style ="none",stub.border.style ="none",row_group.border.top.style ="solid",row_group.border.top.width =px(1),row_group.border.bottom.style ="none",row.striping.include_table_body =FALSE,table.font.size =px(12),table.font.names ="Times New Roman",heading.title.font.size =px(12),heading.subtitle.font.size =px(12),heading.align ="left",data_row.padding =px(4) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_title("title") ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_title("subtitle") ) %>%tab_style(style =cell_text(align ="left"),locations =cells_stub() ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_row_groups() ) %>%##### bolding significancetab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(M1_b_ci, M1_p), rows = sig_rows_m1) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(M2_b_ci, M2_p), rows = sig_rows_m2) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(M3_b_ci, M3_p), rows = sig_rows_m3) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(M4_b_ci, M4_p), rows = sig_rows_m4) )}#### tblstbl_com <-build_table_rq2(com_m1, com_m2, com_m3, com_m4_ml,"Communicative Parenting", table_number =1)tbl_ctrl <-build_table_rq2(ctrl_m1, ctrl_m2, ctrl_m3, ctrl_m4_ml,"Controlling Parenting", table_number =2)tbl_hos <-build_table_rq2(hos_m1, hos_m2, hos_m3, hos_m4_ml,"Hostile/Rejection Parenting", table_number =3)```### Results#### Warm/Communicative ParentingModel 1 included individual-level predictors: youth total problems, child age, child sex, and maternal education. SDQ total problems significantly predicted communicative parenting (β = -0.174, SE = 0.079, p = .028), such that youth with higher levels of total problems perceived themselves as receiving less warm and communicative parenting. Child age also emerged as a significant predictor (β = -0.156, SE = 0.079, p = .048), with older youth reporting lower levels of communicative parenting. Child sex and maternal education were not significant predictors of communicative parenting.Model 2 added maternal adversity indicators (maternal ACEs, informal settlement index, and DSM cross-cutting symptoms), none of which significantly predicted communicative parenting (all ps \> .30). The effects of SDQ total (β = -0.184, SE = 0.081, p = .023) remained significant, while child age became marginally non-significant (β = -0.146, SE = 0.079, p = .063).Model 3 added mother social connections, which was not a significant predictor (β = 0.024, SE = 0.084, p = .773). Model 4 added community-level predictors (community participation and perceived mutual aid), neither of which significantly predicted communicative parenting (all ps \> .55). In the final model, only SDQ total problems remained a significant predictor (β = -0.188, SE = 0.081, p = .021), with child age remaining marginally non-significant (β = -0.149, SE = 0.079, p = .060).None of the mother-level or community-level predictors reached significance across any model. These findings are summarized in @tbl-com.```{r}#| label: tbl-com#| tbl-cap: "Pooled Multilevel Model Estimates for Warm/Communicative Parenting"tbl_com```#### Controlling ParentingModel 1 included youth total problems, child age, child sex, and maternal education as predictors. SDQ total problems significantly predicted controlling parenting (β = 0.324, SE = 0.069, p \< .001), such that youth with higher total problems reported higher levels of controlling parenting. Child age (β = -0.213, SE = 0.067, p = .002) and child sex (β = -0.564, SE = 0.129, p \< .001) were also significant, with older youth and female youth reporting lower levels of controlling parenting. Maternal education was not significant (β = -0.080, SE = 0.081, p = .321).Model 2 added maternal adversity indicators, none of which significantly predicted controlling parenting (all ps \> .23). The effects of SDQ total, child age, and child sex remained stable and significant across Models 2 and 3 with minimal change in magnitude. Model 3 added mother social connections, which was not a significant predictor (β = 0.027, SE = 0.086, p = .756). Model 4 added community-level predictors, neither of which significantly predicted controlling parenting (all ps \> .50).In the final model, SDQ total (β = 0.323, SE = 0.070, p \< .001), child age (β = -0.214, SE = 0.067, p = .001), and child sex (β = -0.567, SE = 0.130, p \< .001) remained significant predictors, with no mother-level or community-level predictors reaching significance across any model. These findings are summarized in @tbl-con.```{r}#| label: tbl-con#| tbl-cap: "Pooled Multilevel Model Estimates for Controlling Parenting"tbl_ctrl```#### Hostile/Rejecting ParentingModel 1 included youth total problems, child age, child sex, and maternal education as predictors. SDQ total problems was a significant predictor of hostile/rejecting parenting (β = 0.425, SE = 0.071, p \< .001), such that youth with higher total problems reported higher levels of hostile and rejecting parenting. Child sex showed a marginal trend in Model 1 (β = -0.269, SE = 0.143, p = .059), suggesting female youth may report lower hostile/rejecting parenting. Child age and maternal education were not significant predictors.Model 2 added maternal adversity indicators, none of which significantly predicted hostile/rejecting parenting (all ps \> .30). Child sex reached significance in Model 2 (β = -0.287, SE = 0.144, p = .047) and remained significant in Model 3. SDQ total remained stable and significant across all models. Model 3 added mother social connections, which was not a significant predictor (β = -0.011, SE = 0.078, p = .892). Model 4 added community-level predictors, with community participation showing a near-significant trend in the opposite direction than hypothesized (β = 0.107, SE = 0.073, p = .140) and perceived mutual aid not significant (β = 0.029, SE = 0.078, p = .704).In the final model, SDQ total (β = 0.427, SE = 0.073, p \< .001) and child sex (β = -0.309, SE = 0.144, p = .032) were significant predictors, with no mother-level or community-level predictors reaching significance.```{r}#| label: tbl-hos#| tbl-cap: "Pooled Multilevel Model Estimates for Hostile/Rejecting Parenting"tbl_hos```#### Overall FindingsAcross all three parenting outcomes none of the mother-level or community-level predictors reached statistical significance across any outcome or model, providing no support for Hypotheses 4 through 6 regarding the role of individual maternal characteristics, relational resources, or community factors in predicting parenting quality. These null findings should be interpreted with caution given the study was underpowered to detect small effects at Level 2, with only 100 mothers at the cluster level.It is also worth noting the ICC pattern across outcomes — controlling parenting showed substantially higher between-family clustering (ICC = 0.355–0.382) compared to communicative and hostile/rejecting parenting (ICC \< 0.01), suggesting that controlling parenting behavior is more consistent within families across siblings than the other parenting dimensions. In other words, warm/communicative and hostile/rejecting parenting seems to be more influenced by youth-level characteristics whereas controlling parenting seems to be driven by family-level factors, albeit none of those included in this study being significant. This might be a function of the underpowered nature of this study, measurement error, or the need to identify other family-level factors that could explain between family variance in controlling parenting.Sequential model comparisons using the D3 pooled likelihood ratio test indicated that no block of predictors significantly improved model fit beyond the previous model for any parenting outcome. These results indicate that neither individual maternal characteristics, relational resources, and community-level factors explained significant additional variance in parenting quality beyond youth characteristics, again, providing no support for Hypotheses 4 through 7.### Sensitivity AnalysesI used 500 iterations of Monte Carlo simulations to randomly select one youth per family. This resulted in a set of 500 subsamples with independent observations.```{r}#|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 in1: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 @fig-resampling.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 @tbl-resampling 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, clustering assumptions, or normality violations observed at Level 1.Given that the study was underpowered to detect small effects at Level 2 due to the sample size of 100 mothers, non-significant results should be viewed with caution. The OLS subsampling analyses corroborate the directional consistency of the estimates; however, the substantial standard deviations observed for several predictors highlight the expected instability in the estimates associated with a limited Level 2 sample.```{r}#| label: fig-resampling#| fig-width: 10#| fig-height: 10#| echo: false#| include: true#| fig-cap: "Distribution of OLS Coefficients Across 500 Monte Carlo Subsamples"all_coefs %>%filter(predictor !="(Intercept)") %>%# remove intercept for claritymutate(outcome =factor(outcome, levels =c("Warm/Communicative", "Controlling", "Hostile/Rejecting")),predictor = dplyr::recode(predictor, !!!predictor_labels_rq2)) %>%ggplot(aes(x = estimate, fill = outcome)) +geom_density(alpha =0.4) +facet_wrap(~ predictor, scales ="free") +scale_fill_manual(values =c("Warm/Communicative"="#F8766D","Controlling"="#00BA38","Hostile/Rejecting"="#619CFF" )) +geom_vline(xintercept =0, linetype ="dashed", color ="black") +labs(x ="Coefficient Estimate", y ="Density", fill ="Outcome") +theme_minimal()``````{r}#| label: tbl-resampling#| tbl-cap: "Sensitivity Analysis: OLS Monte Carlo Subsampling vs. MLM Pooled Estimates"#| echo: falsepredictor_labels <-c("sdq_total"="SDQ Total","child_demo_age"="Child Age","child_demo_gender1"="Child Sex","adult_education"="Maternal Education","mother_aces"="Maternal ACEs","informal_settlement_index"="Informal Settlement Index","mother_dsmtotal"="DSM Cross-Cutting Symptom","mother_sc"="Mother Social Connections","commparticipation"="Community Participation","mutual_aid"="Perceived Mutual Aid")# sim summarysim_summary <- all_coefs %>%filter(predictor !="(Intercept)") %>%mutate(predictor = dplyr::recode(predictor, !!!predictor_labels)) %>%group_by(outcome, predictor) %>%summarise(ols_mean =sprintf("%.3f", mean(estimate)),ols_sd =sprintf("%.3f", sd(estimate)),.groups ="drop" )# extracting pooled estimaetes for both approachesextract_mlm_for_table <-function(model, outcome_label) { est <-testEstimates(model, extra.pars =TRUE)$estimates est <-as.data.frame(est) est$term <-rownames(est) pval_col <-grep("^P", names(est), value =TRUE)[1] est <- est[!is.na(est[[pval_col]]), ] est <- est[est$term !="(Intercept)", ] est$mlm_b <-sprintf("%.3f", est$Estimate) est$mlm_se <-sprintf("%.3f", est$Std.Error) est$mlm_p <-ifelse(est[[pval_col]] < .001, "< .001", sprintf("%.3f", est[[pval_col]])) est$sig <- est[[pval_col]] < .05 est %>%mutate(predictor = dplyr::recode(term, !!!predictor_labels),outcome = outcome_label ) %>%select(outcome, predictor, mlm_b, mlm_se, mlm_p, sig)}mlm_summary <-bind_rows(extract_mlm_for_table(com_m4_reml, "Warm/Communicative"),extract_mlm_for_table(ctrl_m4_reml, "Controlling"),extract_mlm_for_table(hos_m4_reml, "Hostile/Rejecting"))# combine both result typescombined <-left_join(sim_summary, mlm_summary, by =c("outcome", "predictor")) %>%mutate(outcome =factor(outcome, levels =c("Warm/Communicative", "Controlling", "Hostile/Rejecting"))) %>%arrange(outcome, predictor)## significancesig_rows <-which(combined$sig ==TRUE)combined <- combined %>%select(-sig)#### gt tablecombined %>%gt(rowname_col ="predictor", groupname_col ="outcome") %>%tab_spanner(label ="OLS Subsampling (n = 500)", id ="ols_span", columns =c(ols_mean, ols_sd)) %>%tab_spanner(label ="MLM Pooled Estimate", id ="mlm_span", columns =c(mlm_b, mlm_se, mlm_p)) %>%cols_label(ols_mean ="Mean",ols_sd ="SD",mlm_b ="β",mlm_se ="SE",mlm_p ="p" ) %>%tab_stubhead(label ="Predictor") %>%tab_source_note(source_note ="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)." ) %>%tab_options(table.border.top.style ="solid",table.border.top.width =px(2),table.border.bottom.style ="solid",table.border.bottom.width =px(2),heading.border.bottom.style ="solid",heading.border.bottom.width =px(1),column_labels.border.top.style ="none",column_labels.border.bottom.style ="solid",column_labels.border.bottom.width =px(1),table_body.border.top.style ="none",table_body.border.bottom.style ="none",stub.border.style ="none",row_group.border.top.style ="solid",row_group.border.top.width =px(1),row_group.border.bottom.style ="none",row.striping.include_table_body =FALSE,table.font.size =px(12),table.font.names ="Times New Roman",heading.title.font.size =px(12),heading.subtitle.font.size =px(12),heading.align ="left",data_row.padding =px(4) ) %>%tab_style(style =cell_text(weight ="bold"),locations =cells_title("title") ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_title("subtitle") ) %>%tab_style(style =cell_text(align ="left"),locations =cells_stub() ) %>%tab_style(style =cell_text(style ="italic"),locations =cells_row_groups() ) %>%# Bold significant MLM estimates onlytab_style(style =cell_text(weight ="bold"),locations =cells_body(columns =c(mlm_b, mlm_p), rows = sig_rows) )```
