---
title: "Well being statistics"
shorttitle: "Well Stats"
author:
- name: Josef Mana
orcid: 0000-0002-7817-3978
corresponding: true
affiliations:
- name: Institute of Psychology, Czech Academy of Sciences
department: Department of Cognitive Psychology
address: Pod Vodárenskou věží 4
city: Prague
country: Czech Republic
postal-code: 18200
bibliography: references.bib
format:
apaquarto-html:
css: muni-style.css
apaquarto-docx: default
self-contained: true
warning: false
echo: false
---
# Methods
## Statistical analyses
All demographic and outcome variables were described via contingency tables for nominal (gender and education level), and ordinal (single item responses) variables and via sample mean ± standard deviation for continuous variables (age, years of teaching experience and questionnaires' sum scores). All descriptive statistics were computed separately for each city size. Furthermore, null distributions of nominal and ordinal variables across city sizes were tested via Pearson's $\chi^2$ test of the null hypothesis that the joint distribution of the cell counts is the product of the row and column marginals. The null hypotheses of zero difference between means of continuous variables across city sizes as well as across education levels were tested via one-way ANOVAs with Type I sum of squares. The null hypotheses of zero difference between means of continuous variables in men and women were tested via two sample t-test with Welch approximation of degrees of freedom and two-sided alternative hypothesis. Age was kept continuous rather than binned; its association with each psychological outcome was therefore quantified by Pearson's correlation coefficient and by the slope of a simple linear regression of the outcome on age, each with its 95% CI.
To establish internal consistency of methods used, Cronbach's coefficient $\alpha$ [@cronbach1951] was computed together with its 95% confidence interval (CI) via @feldt1987 procedure. Methods were considered to show adequate internal consistency if the $\alpha$ exceeded 0.80.
Reliability of each scale was characterised further by the standardised Cronbach's $\alpha$, by Guttman's $\lambda_6$, i.e., the lower bound based on the squared multiple correlations, and by the average inter-item correlation, all computed after reverse-coded items had been rescored. Because all items are ordinal, and because $\alpha$ additionally assumes tau-equivalence, confirmatory factor analyses (CFA) were fitted with the items treated as ordered categorical, estimated by the diagonally weighted least squares estimator with robust standard errors and a mean- and variance-adjusted test statistic (WLSMV) and with the residuals parametrised, via the lavaan [@rosseel2012] R package. Each questionnaire was fitted both as the unidimensional model implied by using its total sum score and under the factor structure published for that instrument: two subscales, school connectedness and teaching efficacy, for the Teacher Subjective Wellbeing Questionnaire, and nine four-item facets for the Job Satisfaction Survey, the latter fitted as nine correlated factors and as nine factors loading on a single second-order factor. Each subscale was additionally fitted on its own. The two stigma scales are published as unidimensional and were fitted as such. Ordinal omega was computed from the completely standardised solution of each model as $\omega = (\sum\lambda)^2 / ((\sum\lambda)^2 + \sum\theta)$, where $\lambda$ are the standardised loadings and $\theta$ the standardised residual variances, and is reported per first-order factor. Model fit was judged by the scaled $\chi^2$ and by the robust CFI, TLI and RMSEA together with the SRMR.
Inspired by the model proposed by @vogel2007, we conducted a moderated mediation analysis of the effect of Well being and Job Satisfaction on Self Stigma either directly or via the Public Stigma. Specifically, we assumed that Well being and Job Satisfaction are causes of Public Stigma perception which in turn is a cause of Self Stigma perception. Moreover, we assume Well being and Job Satisfaction have a direct causal effect on Self Stigma perception and that all these relations are moderated by age and gender. Finally, Well being and Job Satisfaction were assumed to have common unobserved causes. These causal assumptions are depicted in @fig-dag in the form of a directed acyclic graph (DAG). To estimate the direct and indirect effects implied by this DAG, a structural equation model (SEM) following @fig-dag was fitted and evaluated via the lavaan [@rosseel2012] R package. Full model specification is available at [https://github.com/josefmana/StigHelp.git](https://github.com/josefmana/StigHelp.git).
{#fig-dag}
Because the assumption that Well being and Job Satisfaction are causes of the two stigmas is open to dispute, the data were re-analysed under two alternative sets of causal assumptions. In the first alternative (@fig-dag-reverse), the direction of these relations is reversed: Public Stigma is assumed to be a cause of Self Stigma as well as of Job Satisfaction and Well being, and Self Stigma is assumed to mediate part of the effect of Public Stigma on the latter two. In the second alternative (@fig-dag-confounded), no causal relation between the stigmas on the one hand and Job Satisfaction and Well being on the other hand is assumed at all, and their association is attributed entirely to unmeasured common causes, represented by residual covariances. Both alternatives retain the core assumption that Public Stigma is a cause of Self Stigma, keep age and gender as moderators of all structural paths, and keep the assumption that Job Satisfaction and Well being share unmeasured common causes. Note that under either alternative, Job Satisfaction and Well being no longer lie on a causal path between the two stigmas, so they are not adjusted for when the effect of Public Stigma on Self Stigma is estimated.
{#fig-dag-reverse}
{#fig-dag-confounded}
The structural cores of the three models, i.e., the models with the moderators left out, are all just-identified over the same four variables. They are therefore observationally equivalent. In other words, they imply exactly the same model-implied covariance matrix and consequently the same log-likelihood and the same information criteria (@tbl-modcomp). These data cannot decide between the three sets of causal assumptions and the two alternatives are accordingly reported as a sensitivity analysis of the primary model.
Finally, correlation structure of questionnaires on levels of sum scores and single items was explored. For item-level analysis, all reverse-coded items were rescored such that higher scores of all items within a questionnaire implied higher level of construct being measured, and then pairwise correlation of all item pairs across all questionnaires were estimated via polychoric correlation coefficient [@olsson1979]. Sum scores' correlations were estimated via Pearson's correlation coefficient and null hypotheses of *r* = 0 were tested via t-test based on formula
$t = \frac{r \sqrt(n-2)}{\sqrt(1-r^2)}$
with standard error:
$se = \sqrt(\frac{1-r^2}{n-2})$,
where *r* = Pearson's correlation coefficient and *n* = sample size.
Across analyses, missing values were addressed via complete cases analysis and tests with p \< .05 were considered "statically significant." Analyses were conducted in language for statistical computing R (version `r with(version, paste(major, minor, sep = '.'))`) via "psych" package [@psych], and correlation matrices were visualised via "corrplot" package [@corrplot].
```{r}
#| label: import
#| output: false
library(targets)
library(tidyverse)
library(gt)
upstore <- here::here("_targets")
tar_source()
tar_source("../../R")
new_folder("tables")
# Read outcome data for the report:
n1 <- tar_read(N, store = upstore)$n1 |> c() # Numbers based on city size
n2 <- tar_read(N, store = upstore)$n2 |> c() # Numbers based on gender
n3 <- tar_read(N, store = upstore)$n3 |> c() # Numbers based on education level
tab1 <- tar_read(main_description, store = upstore)
tab2 <- tar_read(education_comparisons, store = upstore)$stats
tab3 <- tar_read(mediation_table, store = upstore)
tab4 <- tar_read(main_correlations, store = upstore)
gtab <- tar_read(gender_comparisons, store = upstore)$stats
alph <- tar_read(alphas_table, store = upstore)
```
# Results
## Sample description
In total, the sample included `r sum(n1)` participants out of which `r n1[1]` came from a city with less than 3,000 inhabitants, `r n1[2]` came from a city with 3,000-15,000 inhabitants, `r n1[3]` came from a city with 15,000-50,000 inhabitants, `r n1[4]` came from a city with 50,000-100,000 inhabitants, `r n1[5]` came from a city with 100,000-500,000 inhabitants, and `r n1[6]` came from a city with more than 500,000 inhabitants. No statistically significant difference in distribution of nominal variables or means of continuous variables between city sizes was observed (@tbl-desc).
```{r}
#| label: tbl-desc
#| tbl-cap: "Sample description"
tab1 <- tab1 |>
mutate(Variable = case_when(
startsWith(Variable, "7.") ~ "Gender",
startsWith(Variable, "9.") ~ "Education level",
startsWith(Variable, "6.") ~ "Age (years)",
startsWith(Variable, "10.") ~ "Teaching experience (years)",
.default = Variable
)) |>
add_column(type = c(rep("Demographic variables", 4), rep("Psychological outcomes", 4))) |>
gt(rowname_col = "Variable", groupname_col = "type") |>
tab_spanner(label = "Number of inhabitants", columns = ends_with("000")) |>
cols_label(stats = html("Statistical analysis")) |>
tab_footnote(
footnote = "presented as women/men/unanswered",
locations = cells_stub(rows = "Gender")
) |>
tab_footnote(
footnote = "presented as High school/Bachelors/Masters/PhD",
locations = cells_stub(rows = "Education level")
) |>
tab_source_note("All continuous variables are reported as sample mean ± standard deviation.")
gtsave(tab1, here::here("tables", "tab1_sample_description.docx"))
tab1
```
There were `r n3[1]` participants with High School, `r n3[2]` participants with Bachelor's degree, `r n3[3]` participants with Master's degree, and `r n3[4]` participants with PhD. We did not observe any statistically significant differences in means of psychological outcomes between different education levels (@tbl-edu).
```{r}
#| label: tbl-edu
#| tbl-cap: "Psychological outcomes' means across education level"
tab2 <- tab2 |>
gt(rowname_col = "scale") |>
tab_spanner(label = "Education level", columns = !contains("stats")) |>
cols_label(stats = "One-way ANOVA") |>
tab_source_note(
"BA = Bachelor's degree, MA = Master's degree. All variables are reported as sample mean ± standard deviation"
)
gtsave(tab2, here::here("tables", "tab2_outcomes_description.docx"))
tab2
```
There were `r n2["female"]` women, `r n2["male"]` men, and `r n2["other"]` participants who did not report their gender. We observed a statistically significant difference between mean Self Stigma ($M_{women}$ = `r with(gtab, mean_female[grepl("Self-stigma", scale)])`, $M_{men}$ = `r with(gtab, mean_male[grepl("Self-stigma", scale)])`, t(`r with(gtab, df[grepl("Self-stigma", scale)])`) = `r with(gtab, t[grepl("Self-stigma", scale)])`, p `r with(gtab, p[grepl("Self-stigma", scale)])`) and mean Public stigma ($M_{women}$ = `r with(gtab, mean_female[grepl("others", scale)])`, $M_{men}$ = `r with(gtab, mean_male[grepl("others", scale)])`, t(`r with(gtab, df[grepl("others", scale)])`) = `r with(gtab, t[grepl("others", scale)])`, p = `r with(gtab, p[grepl("others", scale)])`). In both cases, men reported higher mean stigmatisation compared to women. We did not observe a significant difference between men and women in either mean Well-being ($M_{women}$ = `r with(gtab, mean_female[grepl("being", scale)])`, $M_{men}$ = `r with(gtab, mean_male[grepl("being", scale)])`, t(`r with(gtab, df[grepl("being", scale)])`) = `r with(gtab, t[grepl("being", scale)])`, p = `r with(gtab, p[grepl("being", scale)])`) or mean Job satisfaction ($M_{women}$ = `r with(gtab, mean_female[grepl("Job", scale)])`, $M_{men}$ = `r with(gtab, mean_male[grepl("Job", scale)])`, t(`r with(gtab, df[grepl("Job", scale)])`) = `r with(gtab, t[grepl("Job", scale)])`, p = `r with(gtab, p[grepl("Job", scale)])`).
## Age
Associations of the psychological outcomes with age are presented in @tbl-age.
```{r}
#| label: tbl-age
#| tbl-cap: "Association of psychological outcomes with age"
tab10 <- tar_read(age_comparisons, store = upstore)$stats |>
gt(rowname_col = "scale") |>
tab_spanner(label = "Pearson's correlation", columns = c(r, r_95CI)) |>
tab_spanner(label = "Linear regression", columns = c(b, b_95CI, t, df, p)) |>
cols_label(
r = html("r"),
r_95CI = "95% CI",
b = html("b"),
b_95CI = "95% CI",
t = html("t"),
df = "df",
p = html("p")
) |>
tab_source_note(
"Pearson's correlation coefficients and unstandardised slopes of a simple linear regression of
each outcome's raw sum score on age in years, with 95% confidence intervals. The t-test and its
p-value refer to the null hypothesis of a zero slope."
)
gtsave(tab10, here::here("tables", "tab10_age_associations.docx"))
tab10
```
## Internal consistency
All scales showed adequate internal consistency in the current sample according to the Cronbach's $\alpha$ estimates ($\alpha_D$ = `r with(alph, alpha[grepl("D scale", scale)])`, $\alpha_U$ = `r with(alph, alpha[grepl("U scale", scale)])`, $\alpha_S$ = `r with(alph, alpha[grepl('S scale',scale)])`, $\alpha_K$ = `r with(alph, alpha[grepl("K scale", scale)])`) for Well being, Job satisfaction, Self-stigma and Public stigma respectively.
A fuller picture of the reliability of each scale, adding the standardised $\alpha$, Guttman's $\lambda_6$, ordinal $\omega$, the average inter-item correlation and the fit of the competing confirmatory models, is presented in @tbl-reliability. The unidimensional model fitted the two stigma scales adequately but fitted neither the Teacher Subjective Wellbeing Questionnaire nor the Job Satisfaction Survey, for both of which the published multidimensional structure fitted substantially better. The two Teacher Subjective Wellbeing Questionnaire subscales correlated only moderately with one another, which qualifies the interpretation of its total sum score as a single dimension. The nine-facet models of the Job Satisfaction Survey improved fit markedly over the unidimensional model but returned an inadmissible solution, indicating that the nine facets are not all separable in this sample.
```{r}
#| label: tbl-reliability
#| tbl-cap: "Reliability indexes and competing confirmatory factor models of all scales"
tab11 <- tar_read(reliability_table, store = upstore)
gtsave(tab11, here::here("tables", "tab11_reliability.docx"))
tab11
```
## Moderated mediation
Effect estimates from the SEM are presented in @tbl-sem. We observed a statistically significant effect of Public Stigma on Self Stigma across age groups and genders. The model implies a significant effect of Job Satisfaction and Well being on Public Stigma in Women. The estimates of total effect of Job Satisfaction on Self-Stigma were compatible with a null model (i.e., no causal effect). On the other hand, the model implies a statistically significant total effect of Well being on Self-stigma in women, especially in older age groups. This effect was driven primarily by the direct effect instead of the indirect effect via Public Stigma.
```{r}
#| label: tbl-sem
#| tbl-cap: "Effects derived from the structural equation model of observed variables"
gtsave(tab3, here::here("tables", "tab3_moderated_mediation.docx"))
tab3
```
## Alternative causal assumptions
As shown in @tbl-modcomp, the structural cores of the primary model and of the two alternatives fit the data identically, so the choice between them cannot be made on statistical grounds.
```{r}
#| label: import-alternatives
#| output: false
tab5 <- tar_read(model_comparison, store = upstore)
tab6 <- tar_read(reverse_table, store = upstore)
tab7 <- tar_read(reverse_correlations_table, store = upstore)
tab8 <- tar_read(confounded_table, store = upstore)
tab9 <- tar_read(confounded_correlations_table, store = upstore)
```
```{r}
#| label: tbl-modcomp
#| tbl-cap: "Fit of the structural cores of the three sets of causal assumptions"
gtsave(tab5, here::here("tables", "tab5_model_comparison.docx"))
tab5
```
Under the first alternative, in which the stigmas are causes of Job Satisfaction and Well being (@fig-dag-reverse), effect estimates are presented in @tbl-sem-reverse and the residual correlation of the two outcomes in @tbl-cor-reverse.
```{r}
#| label: tbl-sem-reverse
#| tbl-cap: "Effects derived from the structural equation model assuming that the stigmas cause Job Satisfaction and Well being"
gtsave(tab6, here::here("tables", "tab6_reverse_mediation.docx"))
tab6
```
```{r}
#| label: tbl-cor-reverse
#| tbl-cap: "Residual correlations implied by the model assuming that the stigmas cause Job Satisfaction and Well being"
gtsave(tab7, here::here("tables", "tab7_reverse_correlations.docx"))
tab7
```
Under the second alternative, in which the association between the stigmas and the two outcomes is attributed to unmeasured common causes only (@fig-dag-confounded), the single causal effect the model implies is presented in @tbl-sem-confounded and the residual correlations standing in for the unmeasured common causes in @tbl-cor-confounded.
```{r}
#| label: tbl-sem-confounded
#| tbl-cap: "Conditional effect of Public Stigma on Self Stigma implied by the model assuming unmeasured common causes"
gtsave(tab8, here::here("tables", "tab8_confounded_effects.docx"))
tab8
```
```{r}
#| label: tbl-cor-confounded
#| tbl-cap: "Residual correlations standing in for the unmeasured common causes"
gtsave(tab9, here::here("tables", "tab9_confounded_correlations.docx"))
tab9
```
## Correlation structure
Correlation matrix of sum scores is presented in @fig-corsum (see also @tbl-corsum for numeric representation). Both Self-stigma and Stigma due to others correlated statistically significantly with Well-being and Job satisfaction but not with age and years of teaching experience. Interestingly, the correlation seems to be confined to several stigma due to others scale items on the item-level (see @fig-corits)
{#fig-corsum}
# Appendix
```{r}
#| label: tbl-corsum
#| tbl-cap: "Pearson's correlations between primary outcome variables of the study"
tab4 <- tab4 |>
rownames_to_column("x") |>
gt(rowname_col = "x") |>
sub_missing(missing_text = "")
gtsave(tab4, here::here("tables", "tab4_correlation_matrix.docx"))
tab4
```
{#fig-corits}
# References