We used the Involved Grandparenting and Child Well-Being 2007 survey, recruited by GfK National Opinion Polls, which is a nationally representative sample of English and Welsh adolescents aged 11–1631,32,33. The sample resulted from the distribution of surveys to schools. From the 103 randomly selected schools, in which classes were randomly chosen for survey distribution, 70 schools returned the questionnaires (response rate: 68%). Respondents completed the questionnaire in a school classroom, and the original sample included 1566 adolescent31. 89% of the respondents identified as white ethnicity, 3.5% identified as black or afro-Caribbean, and 4.3% identified as Asian or mixed ethnicity (for more details on respondents, please see e.g.18). When filling in the questionnaire on grandparental investment, respondents were asked to answer questions for only those grandparents who were still alive. Hence, only those respondents who had at least one living grandparent (n = 1488) were considered in the analyses. We also excluded those children from the analyses (n = 58) who were co-residing with their grandparents; such cases are unusual in England and Wales, and it is difficult to estimate the level of grandparental investment in three-generational compared to two-generational households. Nineteen grandchildren were further dropped from the analysis because they had no response in any of the questions asked or had missing data in their age (mean = 13.39, s.d. = 1.41). The total number of children included in the analyses was hence 1411. The proportion of living maternal grandmother, maternal grandfather, paternal grandmother, and paternal grandfather were 83.7%, 68.8%, 73.2%, and 57.1%, respectively.
To measure grandparental investment in their grandchildren, we used questions developed by Elder and Conger34. From the list of all questions available, we chose four questions that directly measured grandparental investment; these were: “how often do you see them” (Q15), “their grandparents had looked after them” (Q26), “they could depend on their grandparents” (Q27), and “provided financial assistance or help” (Q38). Question Q26 was reverse-scaled to match the meaning and ordering of the other scales. Questions Q15, Q26, and Q27 were measured on a 4-point Likert-type item ranging from 1 = not at all/never to 4 = a lot/every day, and Q38 was measured on a 3-point Likert-type scale ranging from 1 = never to 3 = usually.
We used Bayesian structural equation modelling (BSEM) with multiple-indicator latent variables35 to simultaneously examine how the survival status (dead or alive) of other grandparents influenced subjects’ investment in their grandchild (Fig. 1). The response variable was an unobserved latent variable representing the construct “grandparental investment” in their grandchild, measured by four effect indicators (i.e., the latent variable was assumed to cause variation in its indicators) that were the questions asked from the grandchildren about their grandparents’ involvement in their life (i.e., Q15, Q26, Q27, and Q38)35. This means that each question contributed differently to how much each grandparent invested in their grandchild. In other words, we did not aim to model each specific component of grandparental investment separately but considered grandparental investment in their grandchildren as an unmeasured (i.e., measured with measurement error) construct including all its components. The question “provided financial assistance or help” was regarded as the most relevant and reliable observed indicator variable of grandparental investment and was thus used as a marker indicator of the latent variable by fixing its unstandardized loading to unity (Fig. 1). All indicators of the latent variable were treated as ordinal with a probit link function. Therefore, the loadings connecting the latent variable to its indicators can be interpreted as the extent to which a one-unit increase in the latent variable score changes the predicted probit index (z-score) in standard deviation units. In SEM with categorical latent variable indicators with probit link, it is assumed that the categories of observed ordinal variables are determined by the thresholds (the number of categories in the observed variable minus one) in the underlying normally distributed latent variable35. These latent variables then become the indicators of the main latent variable (here, grandparent’s investment), which are, in turn, associated with the ordinal observed variables by the respective threshold structure (Fig. 1)35. Note that when the latent variable with discrete indicators is regressed on independent variables, these coefficients are linear regression coefficients. To be able to compare grandparental investment among the grandparent types using latent variables, a certain level of measurement invariance needs to be established between the grandparent types36. For these data, we relied on partial measurement invariance, where one of the four factor loadings was non-invariant between the groups (thresholds were found to be invariant among grandparent types)21. Moreover, as commonly done in dyadic analyses37, we allowed for covariances among the errors (i.e., unexplained parts of the variation) of the latent variables to account for unmeasured factors that influenced grandparental investment within lineages (Fig. 1). Between-lineage error covariances were also modelled because the same grandchild evaluated investment for each grandparent, which may also have induced correlations between the responses (Fig. 1).
We did not have access to direct measures of grandparental wealth or other resources for these data, which are likely important within-lineage confounders of the effect of grandparent’s survival status on his or her spouse’s investment in grandchildren17. The data set used here did include variables like marriage and employment status of grandparents at the time of the research, which could be regarded, to a varying degree, as proxies of within-lineage resources available for grandparents (e.g., money and time). Another likely confounder of the survival status-investment-relationship is grandparental age, as a focal grandparent’s age is likely linked to the survival status of his or her partner (i.e., older grandparents, particularly grandmothers owing to the general survival advantage of women over men, are more likely to be widows) and his or hers investment in grandchildren (i.e., very old grandparents are likely unable to invest much). However, a major drawback of these data is that information on grandparental age, marriage status, and employment status (along with the data on investment in grandchildren) were only recorded for those grandparents who were alive during the study. Including those potential confounders into the model using the default listwise procedure therefore would have invalidated the whole analysis since all the deceased grandparents would have been omitted from the analysis. As grandparent’s survival status can also be regarded as missing data indicators, handling missingness in grandparent-specific covariates by bringing them into the model by estimating their means, variances, and covariance38, or using multiple imputation would have resulted in non-estimable parameters (i.e., covariances among covariates and survival status variables)39. Consequently, we were unable to include grandparent-specific independent variables in the analysis. Only grandchild age was included as a covariate to improve the efficiency of the regression estimates21. The correlation matrix of the variables used in the analysis can be found in the supplementary materials (Tables S1 and S2).
Instead, we performed a robustness check on our base model to evaluate the impact of unobserved confounding (i.e., shared causes for both independent and dependent variables or omitted variable bias) on the association between grandparental investment and the survival status of other grandparent types, as recently described in Harring et al.40 (please also see Imai et al.41). In this method, the effect of a potential unmeasured confounder, or confounders, is mimicked by a phantom variable that affects both the predictor (i.e., survival status of other grandparents) and the outcome (i.e., investment of a focal grandparent). Phantom variables are latent variables without any indicators, precluding the need for actual data. Instead, the mean and variance of the phantom variables are fixed constants, usually set to zero and unity, respectively40. The rationale is to examine the sensitivity of the original conclusions when one adds the phantom variable as an unmeasured confounder(s) into the model and varies the strength of the expected confounding40.
As discussed above, one potential confounder of the within-lineage effect between the focal grandparent’s investment and the survival status of his or her spouse is the socioeconomic status or resource availability of the grandparents. It is likely that high socioeconomic status in grandparents improves both their survival (e.g., via healthier lifestyle) and increases their ability to invest in their grandchildren. On the other hand, grandparental age is likely a within-lineage confounder that has negative effects both on grandparents’ survival and their investment in grandchildren. However, since the consequences of these two confounders (i.e., resource availability having positive and grandparental age having negative effects on both the outcome and the independent variable) on the association of main interest here are quantitatively the same42, we concentrated only on the case of confounding by resource availability in our sensitivity analysis. While the signs of the suspected confounders are usually easy to argue, the strength of these effects is usually arbitrary without strong prior subject knowledge. Two scenarios were evaluated here. First, the level of confounding was assumed to be roughly equal to the maximum within-lineage effect size observed between the survival status of a spouse and focal grandparent’s investment. Second, we doubled this level of confounding.
We applied Bayesian inference using the Gibbs sampler for the Markov Chain Monte Carlo (MCMC) algorithm to draw posterior distribution to our model parameters. The median of posterior distribution was used as a point estimate and the highest posterior density (HPD) was used for 95% (credibility) interval estimation. Missing data in indicators (which were treated as response variables) were assumed to be missing at random and handled by Bayes as a full-information estimator. Non-informative normally distributed priors were used for structural regression coefficients (hyperparameters for prior mean and variances = N(0, 1002)), factor loadings (N(0, 5)), thresholds (N(0, 1)), and non-informative inverse Wishart priors for error variances (IW(1, 5)), and covariances (IW(0, 5)) of latent variables.
Three chains with a total of 300,000 iterations were run, thinned by every 50th iteration due to some strong autocorrelation among the threshold parameters, with a burn-in of 150,000 iterations. The convergence of MCMC chains was determined using a potential scale reduction factor that compared the estimated between-chains and within-chains variances for each parameter43. In general, values below 1.2 and 1.1 are considered to indicate good convergence of the chains. The potential scale reduction factor for our model was 1.002 after the iterations used here, suggesting appropriate convergence. We also inspected the individual trace plots of individual parameters as well as their autocorrelation plots, confirming convergence. Mplus 8.7 was used for all data analyses44.
Ethical approval and consent to participate
This study was approved, and the research was performed in accordance with the guidelines of University of Oxford Research Committee. All the participants and their parents gave a written consent to participate in the study in accordance with the Declaration of Helsinki.