Stage-specific family structure models: implicit parameter restrictions and Bayesian model comparison with an application to prosocial behavior

Abstract

The two most frequently used specifications in stage-specific family structure analyses of young adult outcomes-state × stage and event × stage-impose restrictions on the parameters of the underlying model of child development. The restrictions imposed by the state specification have substantive a priori disadvantages, implying that use of the state specification requires justification by a specification test. Because the state and event specifications are non-nested, a classical approach to specification testing runs into practical problems. We demonstrate the advantages of a Bayesian approach to the specification testing problem using two applications to young adult prosocial behavior-charitable giving and volunteering. Substantive results are that family structure transitions during adolescence are negatively associated with subsequent giving, and transitions during middle childhood are positively associated with subsequent volunteering. There are other substantive results that would have been missed had we not realized that the state and event specifications impose restrictions on the underlying model and had we not made specification testing a prerequisite for substantive analysis.

Appendix: Brief overview of the Bayesian estimation method

For readers less accustomed to Bayesian methods we provide a brief overview, but see Koop (2003) and Geweke (2005) for complete details. In the Bayesian approach the β parameters (our short-hand term for all the parameters in equation (1)) are random variables. Estimation proceeds by selecting a likelihood function for the y (the vector of outcomes y _i , i = 1,…, n), selecting a prior density for the β parameters, and writing down the posterior density of the β parameters using Bayes’ Theorem:

$$ p\left( {{\varvec{\beta}} ,h|{\varvec{y}}} \right) = \frac{{p\left( {{\varvec{y}}|{\varvec{\beta}} ,h} \right)p_{0} \left( {{\varvec{\beta}} ,h} \right)}}{p\left( {\varvec{y}} \right)} $$

(2)

where p(y | β, h) is the likelihood, p ₀(β, h) is the prior density, and the parameter h ≡ 1/σ ² (σ ² is the variance of u _i ). Once the posterior density p(β, h | y) is known, it can be used to estimate any quantity of interest. For example, the mean β—analogous to the point estimate of the slope parameters in ordinary least squares—is:

$$ {\text{E}}\left[ {\varvec{\beta}} \right] \, = \int {{\varvec{\beta}} p\left( {{\varvec{\beta}} ,h|{\varvec{y}}} \right){\text{ d}}{\varvec{\beta}} } $$

(3)

In our models y is either log charitable giving or log volunteer hours. We assume the u _i are normal, hence the likelihood function of the data y is normal conditional on the β parameters and h. For the parameters’ prior density we use the independent normal-gamma prior: p ₀(β, h) = p _β0(β) × p _h0(h) where p _β0(β) is normal and p _h0(h) is Gamma. For most of the β parameters we set the mean of the prior density equal to zero. In a few cases we have good prior information about a parameter and set the prior accordingly (e.g., we set the prior mean elasticity of giving with respect to current young adult income to 1.0). We set the prior mean of h to be $1/\widehat{\sigma }^{2}$ where $\widehat{\sigma }^{2}$ comes from maximum likelihood estimation of the model. We understand that some may prefer different prior means, so accordingly we set the prior variances of the β parameters to be relatively large (about the same magnitude as the prior variance of the outcome $\widehat{\sigma }^{2}$, implying prior standard deviations on the β parameters of around 3.0), and furthermore give the prior density very little weight relative to the data by setting the prior’s degrees of freedom to only five (compared to the data’s n = 1,011).

The posterior density p(β, h | y) cannot be analytically solved, but is simulated using the Gibbs sampler with data augmentation (see Koop 2003 and Geweke 2005; the algorithm was developed by Chib 1992). We take 1,000 burn-in draws from the posterior density to allow the effect of the initial draw to dissipate, and then take 10,000 additional draws to simulate the posterior. We check convergence by making sure that the 10,000 draws are not autocorrelated and by checking Geweke’s convergence diagnostic (Koop 2003, p. 66). Computations were done in part using the Bayesian Analysis, Computation, and Communication software (http://www.econ.umn.edu/~bacc) and in part using our own software. We also check convergence by (i) thinning the chain (getting our 10,000 draws by taking 1 out of every 10 from 100,000 draws) to reduce autocorrelation between draws and (ii) increasing the number of burn-in draws to 50,000.

The starting point for Bayesian model comparison is to re-write equation (2) to show explicit dependence on the particular model specification of family structure being used (M _m ):

$$ p(\left. {{\varvec{\beta}},h} \right|{\varvec{y}},M_{m} ) = \frac{{p(\left. {\varvec{y}} \right|{\varvec{\beta}},h,M_{m} )\,p_{o} (\left. {{\varvec{\beta}},h} \right|M_{m} )}}{{p(\left. {\varvec{y}} \right|M_{m} )}} $$

(4)

where p(y | M _m ) is the “marginal likelihood” conditional on M _m being the correct specification and m = 1,…,5 indexes the candidate family structure specifications. We estimate the five marginal likelihoods using the Gelfand-Dey method (Koop pp. 104–106). We use these estimates and our prior beliefs about the probability that M _m is the correct model ((p(M _m ) = 1/5 for all m) in Bayes Theorem (again) to calculate the posterior odds ratio between any pair of models—e.g., p(M ₂|y)/p(M ₁ | y) – where p(M _m | y) is the posterior model probability that M _m is the correct model. An assumption that our set of models is exhaustive:

$$ \sum\limits_{m = 1}^{5} {p(\left. {M_{m} } \right|{\varvec{y}}) = 1} $$

allows us to calculate the five posterior model probabilities.

In some applications the posterior model probabilities may indicate that one model specification is strongly preferred: p(M _m |y) ≫ p(M _m′ | y), m ≠ m′. In other applications it may be that more than one model has a non-negligible posterior model probability. However, when such ambiguity arises it can be handled in an intuitively appealing way using Bayesian model averaging: the posterior of any quantity of interest is a weighted average (weights = p(M _m | y)) of the posterior densities p(β, h | y, M _m ) from the models with a non-negligible probability.