Doubly robust estimation of multivariate fractional outcome means with multivalued treatments

Abstract.

This article suggests a doubly robust method of estimating potential outcome means for multivariate fractional outcomes when the treatment of interest is unconfounded and can take more than two values. The method involves maximizing a propensity score weighted multinomial quasi-log-likelihood function with a multinomial logit conditional mean. We show that this estimator, which we call weighted multivariate fractional logit (wmflogit), consistently estimates the potential outcome means if either the propensity score model or the conditional mean model is misspecified. Our simulations demonstrate this double robustness property for the case of shares generated using a Dirichlet distribution. Finally, we advocate for the use of wmflogit by applying it to estimate time-use shares of women participating in the Mexican conditional cash transfer program, Progresa, using Stata’s fmlogit command developed by Buis.

Keywords:

• Double robustness
• fractional outcomes
• multinomial
• unconfoundedness

JEL codes::

• C21
• C31

1. Introduction

In applied work, it is common to encounter outcomes measured as shares or proportions from some total such as pension plan participation rates, expenditure shares, or proportion of time allocated to different household activities. This means that rather than observing a scalar outcome, one obtains a vector of shares that are each bounded in the unit interval and sum to one. Because of their restricted nature, researchers frequently use the multivariate fractional logit estimator (mflogit), proposed by Mullahy (2015), for modeling the conditional mean of multivariate fractional outcomes.

When interest lies in causal inference, one drawback of the mflogit estimator is that consistent estimation of treatment effects (ATEs and ATTs), defined using unconditional potential outcome (PO) means, depends crucially on correct multinomial logit specification of the conditional mean. In this article, we relax this requirement by proposing a propensity score (PS) weighted multivariate fractional logit estimator (wmflogit) that is consistent for the PO means if either the multinomial logit mean or the PS model is correctly specified, but does not require both to be correct at the same time. This makes wmflogit doubly robust (DR) for the PO means.

Our framework considers the general case of multivalued treatments with more than two treatment levels. Such a generalization is useful in a variety of contexts, such as estimating returns to years of education (Uysal, 2015), evaluating impacts of different programs bundled together in some over-arching policy (Frölich, 2004), and even experiments (Negi and Wooldridge, 2020). We use this framework to describe the PO means for multivariate share data. An implication of estimating all PO means in a DR manner is that all pairwise treatment effects, corresponding to any two levels g and h, are also DR.

This article contributes to two different strands of the econometrics literature. The first is to multivariate fractional models, which have been employed frequently in modeling multiple share equations. Linear models rarely provide a good approximation to the conditional mean of share data. These neither guarantee predicted values in the desired range nor are able to account for other specialized features, such as observing outcomes on the boundary. In addition, these do not satisfy the adding up restriction with more than two outcome levels. For fractional outcomes that do not take values on the boundary, Papke and Wooldridge (1996) discuss that a common practice is to model the conditional mean of the log-odds ratio as a linear function $$\begin{array}{l}\mathbb{E}\left[\text{log}\left(\frac{y}{1-y}\right)|\mathbf{x}\right]=\mathbf{x}\mathit{\beta }.\end{array}$$

The appeal of this formulation is that the log-odds has unbounded support, making a linear model for its conditional mean plausible. However, the conditional mean of the outcome, y, cannot be recovered from using the log-odds transformation without imposing additional distributional assumptions. Moreover, it requires ad-hoc adjustments when there are boundary values in the data.

In contrast, multivariate fractional logit provides a tractable functional form for the conditional mean of shares by accommodating their boundary and adding-up constraints. Papke and Wooldridge (1996) were the first to propose fractional logit for a univariate response and Mullahy (2015) and Murteira and Ramalho (2016) generalized it to the case of multivariate fractional responses. Estimation of mflogit uses results from the linear exponential family (LEF) of distributions where an appropriate choice of quasi-log-likelihood (QLL) function delivers consistent parameter estimates. The key result assumes correct specification of the mean but allows the density associated with the QLL to be misspecified (Gourieroux et al., 1984).

This article is also related to the broad literature on DR estimators.¹ For example, Cattaneo (2010) and Ao et al. (2021) propose non parametric doubly robust estimators for the ATE and ATT with multivalued treatments. More specifically, our article contributes directly to the literature that studies the DR property using QML with the LEF. This include Wooldridge (2007), Sloczyński and Wooldridge (2018), Uysal (2015), and Negi (2020). For instance, Wooldridge (2007) studies inverse probability weighted M-estimators for missing data problems. One application of his approach is to the study of DR estimators for the average treatment effect which use QML estimation with specific mean functions. Our article is closest to the univariate analysis of Sloczyński and Wooldridge (2018) and Uysal (2015) in the sense that both these papers focus on DR estimators for treatment effect parameters under unconfoundedness using a multivalued treatment framework. The most salient difference between these papers and our article is that we are interested in multivariate fractional outcomes for which weighted linear regression, as advocated in Uysal (2015), is not an attractive option. We recommend using wmflogit which is a more appropriate DR estimator given the restricted nature of the outcomes. It uses another property of quasi-maximum-likelihood (QML) estimation with the LEF which is that if the mean function is chosen to be the canonical link of an associated QLL, then one obtains consistent PO mean estimates despite the link function being misspecified for the conditional mean (see Wooldridge, 2007; Sloczyński and Wooldridge, 2018). In the case of wmflogit, one uses a multinomial logit mean specification embedded in a weighted multinomial QLL to obtain an estimator that is consistent for the PO means.

Our simulations demonstrate this double robustness property quite well for a 4×1 vector of shares generated using the Dirichlet distribution, where we consider three different treatment levels. We compare wmflogit to four other estimators, namely, linear regression adjustment (lra), inverse probability weighting (ipw), weighted linear regression adjustment (ipwlra), and mflogit. We look at four different cases of misspecification that can arise: correct mean with correct weights, misspecified mean with correct weights, correct mean with misspecified weights, and misspecified mean with misspecified weights. We find that wmflogit outperforms all others estimators in terms of root mean squared error (RMSE) in the vast majority of cases. To illustrate the method in empirical settings, we use wmflogit to estimate the average share of time spent on three different time-use activities by women who participated in the Mexican conditional cash transfer program, Progresa.

The rest of the article is as follows. Section 2 discusses the multivariate PO framework with applications to multinomial or fractional multinomial outcomes. Section 3 considers identification of treatment effects under unconfoundedness and discusses the main assumptions. Section 4 introduces the proposed wmflogit estimator and establishes the result on double robustness along with a discussion of its asymptotic properties. Section 5 talks about other estimators such as ipw, lra, and ipwlra as comparisons to wmflogit and why these are not attractive options when dealing with multivariate fractional outcomes. The same section also discusses the special case of a randomized treatment. Section 6 presents the results from a simulation exercise which generates a 4×1 vector of shares using the Dirichlet distribution. Section 7 applies wmflogit to estimating time use shares of women in Progresa, and Section 8 concludes.

2. Multivariate outcomes and multivalued treatments

Let 𝐘(g) be a J×1 vector of potential outcomes corresponding to treatment state $g=1,\dots ,G$ such that $$\begin{array}{l}\sum _{j=1}^J{Y}_j(g)=1\text{ where each }{Y}_j(g)\in [0,1].\end{array}$$

Examples include shares, proportions, or a binary choice from among J mutually exclusive alternatives. Then, one can define the unconditional PO means corresponding to the j-th alternative as (1) $$\begin{array}{l}{\mu }_{jg}=\mathbb{E}[Y_j(g)]\text{, }j=1,\dots ,J\text{ and }g=1,\dots ,G.\end{array}$$(1)

Also, let (2) $$\begin{array}{l}\mathbb{E}[Y_j(g)|W_h=1]\equiv {\mu }_{jg|h}\end{array}$$(2) be the PO mean for treatment state g among those receiving treatment h. We also let $\mathbf{W}=(W_1,\dots ,W_G)$ denote the vector of binary treatment indicators such that $W_1+\dots +W_G=1$. The probability of being assigned to treatment group, g, is positive, $$\begin{array}{l}{\rho }_g\equiv \mathbb{P}(W_g=1)>0\end{array}$$ such that ${\rho }_1+\dots +{\rho }_G=1$. The observed multivariate fractional outcome can then be written as (3) $$\begin{array}{l}\mathbf{Y}=W_1\cdot \mathbf{Y}(1)+\dots +W_G\cdot \mathbf{Y}(G).\end{array}$$(3)

We can then use (1(1) $$\begin{array}{l}{\mu }_{jg}=\mathbb{E}[Y_j(g)]\text{, }j=1,\dots ,J\text{ and }g=1,\dots ,G.\end{array}$$(1) ) and (2(2) $$\begin{array}{l}\mathbb{E}[Y_j(g)|W_h=1]\equiv {\mu }_{jg|h}\end{array}$$(2) ) to define pairwise average treatment effects for any two treatment levels $g,h=1,\dots ,G$ as (4) $$\begin{array}{l}{\tau }_{jgh}={\mu }_{jg}-{\mu }_{jh}\text{ and }{\tau }_{jhg|h}={\mu }_{jh|h}-{\mu }_{jg|h}\end{array}$$(4) where τ_jgh is the alternative-j specific average treatment effect of treatment g relative to treatment h (ATEs) whereas τ_jhg|h is alternative-j specific average treatment effect of g relative to h among the sub population receiving treatment h, also known as ATTs. These alternative-specific treatment effects can be stacked to define the vector of ATEs or ATTs as (5) $$\begin{array}{l}{\mathit{\tau }}_{gh}={\mathit{\mu }}_g-{\mathit{\mu }}_h\text{ and }{\mathit{\tau }}_{hg|h}={\mathit{\mu }}_{h|h}-{\mathit{\mu }}_{g|h}.\end{array}$$(5)

3. Unconfounded treatment

In many studies, a perfectly randomized treatment or intervention is hard to achieve. Let $\mathbf{X}=(X_1,\dots ,X_K)$ be a vector of observed pre-treatment covariates. In this article, we assume that the treatment is random once we condition on these covariates.² Formally,

Assumption 1.

(Unconfounded Assignment). Assignment is independent of the potential outcomes conditional on the observed covariates. $$\begin{array}{l}\mathbf{W}\mathrm{\perp }[Y_j(1),\dots ,Y_j(G)]\mid \mathbf{X},\mathrm{\forall }j=1,2,\dots ,J.\end{array}$$

Define ${\mu }_{jg}(\mathbf{X})=\mathbb{E}[Y_j(g)|\mathbf{X}]$ to be the conditional mean of the g-th PO corresponding to the j-th alternative. Unconfoundedness allows us to recover or identify this conditional potential outcome mean as a function of the observed population vector $(\mathbf{W},{\mathbf{Y}}^{\prime },\mathbf{X})$ such that $$\begin{array}{l}{\mu }_{jg}(\mathbf{X})=\mathbb{E}[Y_j|\mathbf{X},W_g=1].\end{array}$$

Let m_jg(𝐗, 𝜸_g) be a parametric model for $\mathbb{E}[Y_j|\mathbf{X},W_g=1]$. This model is said to be correctly specified for the true mean if for some true parameter value, ${\mathit{\gamma }}_g^0\in {\mathbf{\Gamma }}_g$, $\mathbb{E}[Y_j|\mathbf{X},W_g=1]=m_{jg}(\mathbf{X},{\mathit{\gamma }}_g^0)$.

Assumption 1 also implies that for each $j=1,\dots ,J$ (6) $$\begin{array}{l}\mathbb{P}(W_g=1|Y_j(1),\dots ,Y_j(G),\mathbf{X})=\mathbb{P}(W_g=1|\mathbf{X})\equiv {\rho }_g(\mathbf{X})\end{array}$$(6) such that, ${\rho }_1(\mathbf{X})+\dots +{\rho }_g(\mathbf{X})=1$ where ρ_g(𝐗) is also known as the generalized PS (Imbens, 2000) for treatment level g.

Assumption 2.

(Overlap). ρ_g(𝐗) > 0 for all $g=1,2,\dots ,G$ and for all 𝐱 in the support of 𝐗.

Assumption 2 also rules out ρ_g(𝐗) = 1 since the probabilities have to sum to one. Overlap is needed in addition to unconfoundedness to identify treatment effects. The assumption implies that for all values of 𝐗, the probability of receiving treatment level g is positive. This allows us to identify conditional average treatment effects which can then be averaged over the marginal distribution of 𝐗 to estimate 𝝉_gh. For the case of estimating ATTs (or 𝝉_hg|h), one can get by with a weaker version of overlap. In that case, we only need ρ_h(𝐗) < 1.

Assumption 3.

(Random Sampling). For a non random integer N, (7) $$\begin{array}{l}\{[{\mathbf{W}}_i,{\mathbf{Y}}_i^{\prime }(1),\dots ,{\mathbf{Y}}_i^{\prime }(G),{\mathbf{X}}_i]\text{: }i=1,2,\dots ,N\}\end{array}$$(7) is independent and identically distributed.

The above assumption restricts dependence between draws and assumes that we have access to a random sample from the population of interest. With standard regularity conditions, we can use i.i.d sequences of the random vector above to apply law of large numbers and central limit theorem.

Given unconfoundedness and overlap, consistent estimators of pairwise ATEs and ATTs can be obtained as (8) $$\begin{array}{l}{\hat{\mathit{\tau }}}_{gh}={\hat{\mathit{\mu }}}_g-{\hat{\mathit{\mu }}}_h\text{ and }{\hat{\mathit{\tau }}}_{hg|h}={\hat{\mathit{\mu }}}_{h|h}-{\hat{\mathit{\mu }}}_{g|h},\text{ for any }g\ne h=1,\dots ,G\end{array}$$(8) where the j-th element of ${\hat{\mathit{\mu }}}_g$ and ${\hat{\mathit{\mu }}}_{g|h}$ is obtained as (9) $$\begin{array}{l}{\hat{\mu }}_{jg}=\frac{1}{N}\sum _{i=1}^N{m}_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)\text{ and }{\hat{\mu }}_{jg|h}=\frac{1}{N_h}\sum _{i=1}^N{W}_{ih}\cdot m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g),j=1,\dots ,J\end{array}$$(9) respectively, with ${\hat{\mathit{\gamma }}}_g$ representing a consistent estimator of ${\mathit{\gamma }}_g^0$ and N_h denoting the total number of units receiving treatment level h in the sample.

The next section discusses the proposed weighted multivariate fractional logit estimators of 𝝁_g and 𝝁_g|h that are doubly robust for the PO means. These estimators involve models for both; the propensity scores and the conditional mean of the outcome, and remain consistent if either of the two models are misspecified.

4. Doubly robust weighted multinomial fractional logit

For each treatment level g, we choose the conditional mean function to be (10) $$\begin{array}{l}{m}_{jg}(\mathbf{X},{\mathit{\gamma }}_g)=\frac{\text{exp}[{\alpha }_{jg}+\mathbf{X}{\mathit{\beta }}_{jg}]}{1+\sum _{h=1}^{J-1}\text{exp}[{\alpha }_{hg}+\mathbf{X}{\mathit{\beta }}_{hg}]},j=1,2,\dots ,J\end{array}$$(10) where ${\mathit{\gamma }}_g=({\mathit{\gamma }}_{1g}^{\prime },\dots ,{\mathit{\gamma }}_{Jg}^{\prime }{)}^{\prime }$, ${\mathit{\gamma }}_{jg}=({\alpha }_{jg},{\mathit{\beta }}_{jg}^{\prime }{)}^{\prime }$, and α_Jg = 0 and 𝜷_Jg = 𝟎 for $g=1,\dots ,G$. Since we are dealing with outcomes that are binary or fractional in nature, using a multivariate fractional logit specification for the conditional mean ensures that fitted values lie in the unit interval since $0<m_{jg}(\cdot )<1$ and satisfies the adding-up constraint i.e., $\sum _{j=1}^J{m}_{jg}(\cdot )=1$ for each g. We will call this as the mflogit specification even though it applies to binary outcomes as well.

Also, let p_g(𝐗, 𝜹) be a parametric model for the generalized propensity score for treatment level g which is said to be correctly specified if there exists some true parameter value, 𝜹₀∈Δ, such that ${\rho }_g(\mathbf{X})=p_g(\mathbf{X},{\mathit{\delta }}_0)$.

Then the proposed wmflogit estimator for μ_jg is obtained using (9(9) $$\begin{array}{l}{\hat{\mu }}_{jg}=\frac{1}{N}\sum _{i=1}^N{m}_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)\text{ and }{\hat{\mu }}_{jg|h}=\frac{1}{N_h}\sum _{i=1}^N{W}_{ih}\cdot m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g),j=1,\dots ,J\end{array}$$(9) ) where ${\hat{\mathit{\gamma }}}_g$ solves a propensity score weighted multinomial QLL given by³ (11) $$\begin{array}{l}\sum _{i=1}^N\sum _{j=1}^J\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot Y_{ij}\cdot \text{log}[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)].\end{array}$$(11)

In the equation above, $p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})$ is the estimated PS obtained from solving a multinomial logit problem in a first-step.⁴ Similarly, for the case of μ_jg|h, ${\hat{\mathit{\gamma }}}_g$ will solve (12) $$\begin{array}{l}\sum _{i=1}^N\sum _{j=1}^J\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}{p}_h({\mathbf{X}}_i,\hat{\mathit{\delta }})\cdot Y_{ij}\cdot \text{log}[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)].\end{array}$$(12)

Essentially, wmflogit can be thought of as an inverse probability weighted regression adjustment (ipwra) estimator which uses the multinomial logit conditional mean. Later we will be comparing wmflogit to the ipwra estimator which assumes the conditional means to be linear. To emphasize it’s linear, we term this estimator ipwlra.

We formalize the double robustness result for the stacked vector of wmflogit PO mean estimators, ${\hat{\mathit{\mu }}}_g$ and ${\hat{\mathit{\mu }}}_{g|h}$, in the proposition below.

Proposition 1.

Assume 1, 2, and 3.

1. If the propensity score model is correctly specified i.e. $\hat{\mathit{\delta }}\stackrel{p}{\to }{\mathit{\delta }}_0$ , then $$\begin{array}{ll}{\mathit{\mu }}_g& =\mathbb{E}[{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\text{ and }{\mathit{\mu }}_{g|h}=\mathbb{E}[{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|W_h=1]\end{array}$$

   where ${\mathit{\gamma }}_g^{\ast }$ ’s solve the population version of the QLL given in (4.2) and (4.3), respectively.

2. If the mean model is correctly specified i.e ${\hat{\mathit{\gamma }}}_g\stackrel{p}{\to }{\mathit{\gamma }}_g^0$ , then $$\begin{array}{ll}{\mathit{\mu }}_g& =\mathbb{E}[{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)]\text{ and }{\mathit{\mu }}_{g|h}=\mathbb{E}[{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)|W_h=1]\end{array}$$

   where $\hat{\mathit{\delta }}\stackrel{p}{\to }{\mathit{\delta }}^{\ast }\ne {\mathit{\delta }}_0$ and ${\mathit{\gamma }}_g^0$ ’s solve the population version of the QLL given in (4.2) and (4.3), respectively.

The proof uses the first-order condition corresponding to the intercept in the population QLL and shows how wmflogit identifies 𝝁_g and 𝝁_g|h under both cases; when the propensity score model is misspecified and when the conditional mean model is misspecified.

The first half of double robustness implies that ${\hat{\mathit{\mu }}}_g$ is consistent for 𝝁_g even if we allow for the possibility that ${\hat{\mathit{\gamma }}}_g\stackrel{p}{\to }{\mathit{\gamma }}_g^{\ast }\ne {\mathit{\gamma }}_g^0$ as long as $\hat{\mathit{\delta }}\stackrel{p}{\to }{\mathit{\delta }}_0$. In contrast, the second half allows for the possibility that if $\hat{\mathit{\delta }}\stackrel{p}{\to }{\mathit{\delta }}^{\ast }\ne {\mathit{\delta }}_0$, then ${\hat{\mathit{\mu }}}_g$ is still consistent for PO means as long as ${\hat{\mathit{\gamma }}}_g\stackrel{p}{\to }{\mathit{\gamma }}_g^0$. In other words, when the mflogit model is correctly specified for the true conditional mean, weighting is not needed to guarantee consistent estimation of the unconditional PO means. However, when the mean is misspecified, a correctly specified propensity score model will guarantee consistent estimation of the PO means despite misspecification of 𝐦_g(𝐗, 𝜸_g). The resulting estimator of 𝝁_g is then doubly robust for the corresponding population PO means and consequently, ${\hat{\mathit{\tau }}}_g$ will also be doubly robust for 𝝉_g.

The double robustness result relies on estimation of ${\hat{\mathit{\gamma }}}_g$ using quasi-maximum likelihood in the linear exponential family (LEF) of distributions. We use features of LEF combined with particular mean functions to show that if one chooses a combination of canonical link and quasi log-likelihood functions carefully, then consistent estimation of PO mean is possible despite misspecification in the conditional mean function. As shown in Table 2 of Negi and Wooldridge (2021), this choice depends on the nature of outcomes one wishes to study. For the case of multivariate fractional outcomes, the multinomial-logit canonical link and the multinomial QLL will guarantee this result. The next section uses the theory of two-step estimators to obtain the asymptotic variance of the wmflogit estimator of 𝝁_g and 𝝁_g|h.

4.1. Asymptotic results

Let ${\mathbf{s}}_{i1}^{\ast }\equiv {\mathbf{s}}_{i1}({\mathit{\gamma }}_g^{\ast },{\mathit{\delta }}^{\ast })$ be the (J−1) vector of scores (or first order conditions) corresponding to the wmflogit problem for ${\hat{\mathit{\gamma }}}_g$ which solves (11(11) $$\begin{array}{l}\sum _{i=1}^N\sum _{j=1}^J\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot Y_{ij}\cdot \text{log}[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)].\end{array}$$(11) ) and ${\mathbf{s}}_{i2}^{\ast }\equiv {\mathbf{s}}_{i2}({\mathit{\delta }}^{\ast })$ be the (G−1) vector of scores for the mflogit problem for estimating the propensity score parameter, 𝜹, where each of these scores are evaluated at the probability limits of $({\hat{\mathit{\gamma }}}_g,\hat{\mathit{\delta }})$, respectively. Also, let ${\mathbf{S}}_{1\gamma }^{\ast }=\mathbb{E}[{\mathrm{\nabla }}_{{\gamma }_g}{\mathbf{s}}_{i1}^{\ast }]$, ${\mathbf{S}}_{1\delta }^{\ast }=\mathbb{E}[{\mathrm{\nabla }}_{\delta }{\mathbf{s}}_{i1}^{\ast }]$ and ${\mathbf{S}}_{2\delta }^{\ast }=\mathbb{E}[{\mathrm{\nabla }}_{\delta }{\mathbf{s}}_{i2}^{\ast }]$. Under standard regularity conditions, the influence function representation for $\hat{\mathit{\delta }}$ is (13) $$\begin{array}{l}\sqrt{N}(\hat{\mathit{\delta }}-{\mathit{\delta }}^{\ast })=-({\mathbf{S}}_{2\delta }^{\ast }{)}^{-1}\cdot N^{-1/2}\sum _{i=1}^N{\mathbf{s}}_{i2}^{\ast }+o_p(1)=-N^{-1/2}\sum _{i=1}^N{\mathbf{S}}_{2\delta }^{\ast -1}{\mathbf{s}}_{i2}^{\ast }+o_p(1).\end{array}$$(13)

The influence function (IF) representation for ${\hat{\mathit{\gamma }}}_g$ is obtained from the wmflogit first-order conditions⁵ given by (14) $$\begin{array}{l}\sum _{i=1}^N{\mathbf{s}}_{i1}({\hat{\mathit{\gamma }}}_g,\hat{\mathit{\delta }})=\mathbf{0}.\end{array}$$(14)

Using a mean value expansion around ${\mathit{\gamma }}_g^{\ast }$ and using Eq. (13(13) $$\begin{array}{l}\sqrt{N}(\hat{\mathit{\delta }}-{\mathit{\delta }}^{\ast })=-({\mathbf{S}}_{2\delta }^{\ast }{)}^{-1}\cdot N^{-1/2}\sum _{i=1}^N{\mathbf{s}}_{i2}^{\ast }+o_p(1)=-N^{-1/2}\sum _{i=1}^N{\mathbf{S}}_{2\delta }^{\ast -1}{\mathbf{s}}_{i2}^{\ast }+o_p(1).\end{array}$$(13) ), we obtain (15) $$\begin{array}{ll}\sqrt{N}({\hat{\mathit{\gamma }}}_g-{\mathit{\gamma }}_g^{\ast })& =-({\mathbf{S}}_{1\mathit{\gamma }}^{\ast }{)}^{-1}\cdot N^{-1/2}\sum _{i=1}^N{\mathbf{d}}_i^{\ast }+o_p(1)\end{array}$$(15) where ${\mathbf{d}}_i^{\ast }={\mathbf{s}}_{i1}^{\ast }-{\mathbf{S}}_{1\delta }^{\ast }{\mathbf{S}}_{2\delta }^{\ast -1}{\mathbf{s}}_{i2}^{\ast }$. Then, given that ${\hat{\mathit{\gamma }}}_g$ is both consistent and asymptotically normal, consistency and normality of ${\hat{\mu }}_{jg}$ follows from standard results in econometrics (see Wooldridge, 2010, chapter 12) where $$\begin{array}{ll}\sqrt{N}({\hat{\mu }}_{jg}-{\mu }_{jg})& =\frac{1}{\sqrt{N}}\sum _{i=1}^N[m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)-{\mu }_{jg}]\\ & =\frac{1}{\sqrt{N}}\sum _{i=1}^N[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g^{\ast })-{\mu }_{jg}]+\mathbb{E}[{\mathrm{\nabla }}_{{\gamma }_g}{m}_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g^{\ast })]\cdot \sqrt{N}({\hat{\mathit{\gamma }}}_g-{\mathit{\gamma }}_g^{\ast })+o_p(1)\\ & =\frac{1}{\sqrt{N}}\sum _{i=1}^N\{{\dot{m}}_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g^{\ast })-{\mathbf{M}}_{jg}^{\ast }\cdot {\mathbf{S}}_{1\mathit{\gamma }}^{\ast -1}\cdot {\mathbf{d}}_i^{\ast }\}+o_p(1)\\ & =\frac{1}{\sqrt{N}}\sum _{i=1}^N{\psi }_{jg}({\mathbf{X}}_i)+o_p(1)\text{ where }\mathbb{E}[{\psi }_{jg}({\mathbf{X}}_i)]=0.\end{array}$$

Then stacking all such IF’s for all $j=1,\dots ,J$, $$\begin{array}{l}\sqrt{N}({\hat{\mathit{\mu }}}_g-{\mathit{\mu }}_g)=\frac{1}{\sqrt{N}}\sum _{i=1}^N{\mathit{\psi }}_g({\mathbf{X}}_i)+o_p(1)\text{ where }{\mathit{\psi }}_{\mathit{g}}\mathbf{(}{\mathbf{X}}_{\mathit{i}}\mathbf{)}={\left(\begin{array}{c}{\psi }_{1g}({\mathbf{X}}_i),\dots ,{\psi }_{Jg}({\mathbf{X}}_i)\end{array}\right)}^{\prime }.\end{array}$$

Then, $\text{Avar}[\sqrt{N}({\hat{\mathit{\mu }}}_g-{\mathit{\mu }}_g)]=\mathbb{E}[{\mathit{\psi }}_g({\mathbf{X}}_i){\mathit{\psi }}_g^{\prime }({\mathbf{X}}_i)]$ provided it exists and $\widehat{\text{Avar}}[\sqrt{N}({\hat{\mathit{\mu }}}_g$ $-{\mathit{\mu }}_g)]=\frac{1}{N}\sum _{i=1}^N{\hat{\mathit{\psi }}}_g({\mathbf{X}}_i){\hat{\mathit{\psi }}}_g^{\prime }({\mathbf{X}}_i)$ where ${\hat{\mathit{\psi }}}_g({\mathbf{X}}_i)$ is evaluated at $({\hat{\mathit{\gamma }}}_g$, $\hat{\mathit{\delta }})$.

Let ${\stackrel{\mathbf{~}}{\mathit{s}}}_{i1}^{\ast }$ be the score of the wmflogit problem solved by ${\hat{\mathit{\gamma }}}_g$ in (12(12) $$\begin{array}{l}\sum _{i=1}^N\sum _{j=1}^J\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}{p}_h({\mathbf{X}}_i,\hat{\mathit{\delta }})\cdot Y_{ij}\cdot \text{log}[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)].\end{array}$$(12) ), ${\stackrel{\mathbf{~}}{\mathit{S}}}_{1\gamma }^{\ast }=\mathbb{E}[{\mathrm{\nabla }}_{{\gamma }_g}{\stackrel{\mathbf{~}}{\mathit{s}}}_{i1}^{\ast }]$, and ${\stackrel{\mathbf{~}}{\mathit{S}}}_{1\delta }^{\ast }=\mathbb{E}[{\mathrm{\nabla }}_{\delta }{\stackrel{\mathbf{~}}{\mathit{s}}}_{i1}^{\ast }]$. The influence function representation for ${\hat{\mathit{\gamma }}}_g$ is (16) $$\begin{array}{ll}\sqrt{N}({\hat{\mathit{\gamma }}}_g-{\mathit{\gamma }}_g^{\ast })& =-({\stackrel{\mathbf{~}}{\mathit{S}}}_{1\mathit{\gamma }}^{\ast }{)}^{-1}\cdot N^{-1/2}\sum _{i=1}^N{\stackrel{\mathbf{~}}{\mathit{d}}}_i^{\ast }+o_p(1)\end{array}$$(16) where ${\stackrel{\mathbf{~}}{\mathit{d}}}_i^{\ast }={\stackrel{\mathbf{~}}{\mathit{s}}}_{i1}^{\ast }-{\stackrel{\mathbf{~}}{\mathit{S}}}_{1\delta }{\mathit{S}}_{2\delta }{\mathit{s}}_{i2}^{\ast }$. The IF representation for μ_jg|h is given as $$\begin{array}{ll}\sqrt{N}({\hat{\mu }}_{jg|h}-{\mu }_{jg|h})& =\frac{1}{\sqrt{N}}\sum _{i=1}^N\{\frac{W_{ih}}{{\hat{\rho }}_h}\cdot m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)-{\mu }_{jg|h}\}\\ & =\frac{1}{\sqrt{N}}\sum _{i=1}^N\{\frac{W_{ih}}{{\rho }_h}\cdot m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g^{\ast })-{\mu }_{jg|h}\}+{\mathbf{M}}_{jg|h}^{\ast }\cdot \sqrt{N}({\hat{\mathit{\gamma }}}_g-{\mathit{\gamma }}_g^{\ast })+o_p(1)\\ & =\frac{1}{\sqrt{N}}\sum _{i=1}^N\{\frac{W_{ih}}{{\rho }_h}\cdot m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g^{\ast })-{\mu }_{jg|h}-{\mathbf{M}}_{jg|h}^{\ast }\cdot {\stackrel{\mathbf{~}}{\mathit{S}}}_{1\mathit{\gamma }}^{\ast -1}\cdot {\stackrel{\mathbf{~}}{\mathit{d}}}_i^{\ast }\}+o_p(1)\\ & =\frac{1}{\sqrt{N}}\sum _{i=1}^N{\tilde{\psi }}_{jg}({\mathbf{X}}_i)+o_p(1)\text{ such that }\mathbb{E}[{\tilde{\psi }}_{jg}({\mathbf{X}}_i)]=0\end{array}$$ where ${\mathbf{M}}_{jg|h}^{\ast }=\mathbb{E}[{\mathrm{\nabla }}_{{\gamma }_g}{m}_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g^{\ast })|W_h=1]$. Again, stacking all the IF’s for $j=1,\dots ,J$ we get $$\begin{array}{l}\sqrt{N}({\hat{\mathit{\mu }}}_{g|h}-{\mathit{\mu }}_{g|h})=\frac{1}{\sqrt{N}}\sum _{i=1}^N{\stackrel{\mathbf{~}}{\mathit{\psi }}}_g({\mathbf{X}}_i)+o_p(1).\end{array}$$

Then $\text{Avar}[\sqrt{N}({\hat{\mathit{\mu }}}_{jg|h}-{\mathit{\mu }}_{jg|h})]=\mathbb{E}[{\stackrel{\mathbf{~}}{\mathit{\psi }}}_g({\mathbf{X}}_i){\stackrel{\mathbf{~}}{\mathit{\psi }}}_g^{\prime }({\mathbf{X}}_i)]$ provided it exists and $\widehat{\text{Avar}}[\sqrt{N}({\hat{\mathit{\mu }}}_g-{\mathit{\mu }}_g)]=\frac{1}{N}\sum _{i=1}^N{\widehat{\stackrel{\mathbf{~}}{\mathit{\psi }}}}_g({\mathbf{X}}_i){\widehat{\stackrel{\mathbf{~}}{\mathit{\psi }}}}_g^{\prime }({\mathbf{X}}_i)$ where ${\widehat{\stackrel{\mathbf{~}}{\mathit{\psi }}}}_g({\mathbf{X}}_i)$ is evaluated at $({\hat{\mathit{\gamma }}}_g$, $\hat{\mathit{\delta }})$. Consistency of the asymptotic variance estimator for ${\hat{\mathit{\mu }}}_g$ and ${\hat{\mathit{\mu }}}_{g|h}$ then follows from theorem 4.1 of Newey and McFadden (1994).

5. Other estimators

This section discusses alternative estimators like propensity score weighting, linear regression adjustment, and weighted linear regression adjustment, which can be used for consistently estimating 𝝁_g and 𝝁_g|h. However, as we discuss below, these estimators are not attractive choices when dealing with multivariate fractional outcomes.

5.1. Weighting with propensity score

With unconfoundedness and overlap, a common strategy to estimate PO means is inverse propensity score weighting (ipw). Since ipw estimators do not involve models for the conditional mean, one can obtain these estimators directly using Lemma 3.2 of Sloczyński and Wooldridge (2018) ⁶ where (17) $$\begin{array}{l}{\hat{\mu }}_{jg}=\frac{1}{N}\cdot \sum _{i=1}^N\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot Y_{ij}\text{ and }{\hat{\mu }}_{jg|h}=\frac{1}{N_h}\sum _{i=1}^N\frac{W_{ig}\cdot p_h({\mathbf{X}}_i,\hat{\mathit{\delta }})}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot Y_{ij},\end{array}$$(17) respectively, where $p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})\stackrel{p}{\to }{p}_g(\mathbf{X},{\mathit{\delta }}_0)={\rho }_g(\mathbf{X})$ for all $g=1,\dots ,G$. In other words, for propensity score weighting to produce consistent estimators of the PO means, it is necessary that the propensity score model is correctly specified. This is again restrictive, as it offers practitioners only one possibility of obtaining consistent PO mean estimates. As a practical matter, the estimators in (17) are limited because the weights do not sum to unity, although this is easy to fix (see Sloczyński and Wooldridge, 2018).

5.2. Linear regression adjustment

Another common approach is to use linear regression adjustment (lra) to model the PO means. In this case, one assumes the conditional means to be linear in parameters i.e., (18) $$\begin{array}{l}{m}_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)={\alpha }_{jg}+{\mathbf{X}}_i{\mathit{\beta }}_{jg},j=1,2,\dots ,J.\end{array}$$(18)

It is also useful to express the conditional means as $m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)={\mu }_{jg}+{\dot{\mathit{X}}}_i{\mathit{\beta }}_{jg}$ where ${\dot{\mathit{X}}}_i={\mathbf{X}}_i-\mathbb{E}({\mathbf{X}}_i)$ are the population-demeaned covariates. This representation makes it clear how Eq. (9(9) $$\begin{array}{l}{\hat{\mu }}_{jg}=\frac{1}{N}\sum _{i=1}^N{m}_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)\text{ and }{\hat{\mu }}_{jg|h}=\frac{1}{N_h}\sum _{i=1}^N{W}_{ih}\cdot m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g),j=1,\dots ,J\end{array}$$(9) ) would give us ${\hat{\mu }}_{jg}$ and ${\hat{\mu }}_{jg|h}$ if we estimate regressions $$\begin{array}{l}{Y}_{ij}\text{ on }1,{\ddot{\mathit{X}}}_i\text{ using }{W}_{ig}=1,\text{ for all }j,g\end{array}$$ where ${\ddot{\mathit{X}}}_i={\mathbf{X}}_i-\overline{\mathit{X}}$ are the sample-demeaned covariates. We can then stack these estimates for all j to obtain ${\hat{\mathit{\mu }}}_g$ and ${\hat{\mathit{\mu }}}_{g|h}$ which can then consequently be plugged in Eq. (8(8) $$\begin{array}{l}{\hat{\mathit{\tau }}}_{gh}={\hat{\mathit{\mu }}}_g-{\hat{\mathit{\mu }}}_h\text{ and }{\hat{\mathit{\tau }}}_{hg|h}={\hat{\mathit{\mu }}}_{h|h}-{\hat{\mathit{\mu }}}_{g|h},\text{ for any }g\ne h=1,\dots ,G\end{array}$$(8) ) to get estimates of pairwise ATEs and ATTs.

It is important to note that the lra estimator of PO means will be consistent only if the conditional mean vector of the outcomes is truly linear in parameters. As we mention in the introduction, this is hard to justify for multivariate fractional responses. An issue with linear models is that they are rarely appropriate for fractional outcomes much like how linear probability models are ill-suited for studying binary outcomes. In addition, lra does not satisfy the adding-up constraint when G > 2. More importantly, since it relies on only one feature of the conditional distribution (the mean), it also does not share the double robustness property of wmflogit, thereby making it less desirable in the current setting.

5.3. Linear regression adjustment with weighting

A method that overcomes the shortcomings of PS weighting and linear regression adjustment is weighted linear regression adjustment (ipwlra). Much like lra, this estimator also assumes a linear conditional mean function. However, ipwlra estimators of ${\hat{\mu }}_{jg}$ and ${\hat{\mu }}_{jg|h}$ solve propensity score weighted linear regressions (19) $$\begin{array}{l}\underset{({\mathit{\mu }}_g,{\mathit{\beta }}_g)}{\text{arg min}}\sum _{i=1}^N\left(\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\right)(Y_{ij}-{\mu }_{jg}-{\ddot{\mathit{X}}}_i{\mathit{\beta }}_{jg}{)}^2\end{array}$$(19) and (20) $$\begin{array}{l}\underset{({\mathit{\mu }}_{g|h},{\mathit{\beta }}_{g|h})}{\text{arg min}}\sum _{i=1}^N(W_{ig}\cdot \frac{p_h({\mathbf{X}}_i,\hat{\mathit{\delta }})}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})})(Y_{ij}-{\mu }_{jg|h}-{\stackrel{\mathbf{˘}}{\mathit{X}}}_{ih}{\mathit{\beta }}_{jg|h}{)}^2\end{array}$$(20) respectively, where ${\stackrel{\mathbf{˘}}{\mathit{X}}}_{ih}={\mathbf{X}}_i-{\overline{\mathit{X}}}_h$, ${\overline{\mathit{X}}}_h=N_h^{-1}\sum _{i=1}^N{W}_{ih}\cdot {\mathbf{X}}_i$, and $p_g(\cdot ,\hat{\mathit{\delta }})$, $g=1,2,\dots ,G$ is an estimator of ρ_g(𝐗). In principle, ipwlra is also robust to misspecification of the propensity score or the outcome model, as long as both are not misspecified at the same time. But because a linear functional form for the mean is rarely appropriate for fractional data, ipwlra is most likely going to rely on correct specification of the propensity scores to obtain a consistent linear approximation to the mean. In contrast, it makes more sense to think that wmflogit, which models the conditional mean as a multinomial logit, will provide a better approximation to the true conditional mean function and consequently has a better chance of being truly DR.

5.4. Randomized experiment

When the vector of treatment indicators, 𝐖, is randomized, 𝐖 is independent of $[\mathbf{Y}(1),\dots ,\mathbf{Y}(G),\mathbf{X}]$ – the vectors of potential outcomes and the covariates 𝐗. In this case, the subsample averages for each treatment level are both unbiased and consistent estimators of the PO means, 𝝁_g (which are the same as the PO means conditional on different treatment levels by randomization). As shown by Negi and Wooldridge (2020) when the potential outcomes 𝐘(g) are scalars, linear regression adjustment separately for each treatment preserves consistency – via a linear projection argument – and is generally asymptotically more efficient than the subsample averages if 𝐗 helps to predict the potential outcomes. Here, the potential outcomes 𝐘(g) are vectors, but it follows immediately from Negi and Wooldridge (2020) that linear regression adjustment applied to each outcome Y_j and treatment level g consistently estimates μ_jg and is generally more efficient than using the subsample averages. One unsatisfying feature of applying linear RA to multinomial or multivariate fractional outcomes is that the fitted values are not guaranteed to be in the unit interval and the adding-up restriction – that fitted values should sum to unity – will not hold. By contrast, the multinomial logit (or fractional logit) estimator will produce fitted values that satisfy the logical restrictions. Importantly, our double robustness result implies that the mflogit estimator is consistent under random assignment even if the conditional means are arbitrarily misspecified. The logic is simple: under random assignment, the weighted mflogit estimator is equivalent to the unweighted estimator because the propensity scores do not depend on 𝐗, and so using separate weighted estimation for each g is the same as using separate mflogit estimation for each g.

In the scalar outcome case, Negi and Wooldridge (2020) show that if the conditional mean is correctly specified, then non linear regression adjustment that uses the canonical link function in the linear exponential family is asymptotically more efficient than linear regression adjustment. We will not show formally what happens in the current case where 𝐘(g) is a vector that satisfies adding-up restrictions. Practically, the important point is that, under random assignment, mflogit is very attractive because it is fully robust to conditional mean misspecification, ensures that fitted values are in the unit interval and sum to one, and is likely to be more efficient – perhaps substantially more efficient – than linear regression adjustment (which is in turn more efficient than subsample averages). See Negi and Wooldridge (2020) for simulation evidence using logistic regression and Poisson regression with an exponential mean.

6. Monte Carlo simulations

This section studies the finite-sample bias and standard deviation of wmflogit in comparison to ipw, lra, ipwlra, and mflogit for estimating PO means of a 4×1 vector of shares (J = 4). These are generated through a Dirichlet distribution where the true conditional mean is multinomial with a quadratic index. We consider sample sizes of $N\in \{5000,8000\}$ and three different treatment levels (G = 3).⁷ So, $$\begin{array}{l}{\mathbf{Y}}_i=W_{i1}{\mathbf{Y}}_i(1)+W_{i2}{\mathbf{Y}}_i(2)+W_{i3}{\mathbf{Y}}_i(3)\end{array}$$ where ${\mathbf{Y}}_i(g)={(Y_{i1}(g),Y_{i2}(g),Y_{i3}(g),Y_{i4}(g))}^{\prime }$ is the outcome-vector of shares for treatment level, $g=1,2,3$. The population vector $\{({\mathbf{W}}_i,{\mathbf{Y}}_i^{\prime },{\mathbf{X}}_i):i=1,2,\dots ,N\}$ is generated using a hundred thousand observations from which samples are drawn without replacement. The empirical distributions are then obtained from 1,000 monte carlo draws corresponding to each sample size. We study the behavior of wmflogit with respect to the other estimators under four realistic estimation scenarios: correct mean and weight, misspecified mean with correct weight, correct mean with misspecified weight, and misspecified mean and weight. These are enumerated in Table 1.

Table 1. Estimation scenarios for simulations.

Case	Assumed mean	Assumed PS	lra	ipw	ipwlra	mflogit	wmflogit
1	MNFL, quadratic	MNL, quadratic	inconsistent	consistent	consistent	consistent	consistent
2	MNFL, linear	MNL, quadratic	inconsistent	consistent	consistent	inconsistent	consistent
3	MNFL, quadratic	MNL, linear	inconsistent	inconsistent	inconsistent	consistent	consistent
4	MNFL, linear	MNL, linear	inconsistent	inconsistent	inconsistent	inconsistent	inconsistent

^aThe true conditional mean is MNFL with index quadratic in the covariates. Therefore, the mean is correctly specified in Cases 1 and 3. The mean is misspecified in Cases 2 and 4. Likewise, because the true PS is generated as MNL with index quadratic in the covariates, the PS is correctly specified in Cases 1 and 2 and misspecified in Cases 3 and 4.

^blra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

Given that the true DGP generates share data using a Dirichlet distribution with multinomial mean, it is more appropriate to think that linear regression and weighted linear regression provide the best linear approximation for multivariate fractional outcomes. In that sense, lra will always be inconsistent for the PO mean vector, 𝝁_g, whereas ipwlra can only be consistent if the PS model is correctly specified. The latter argument also holds for ipw which explicitly depends on correct weights for identification of the PO means given an unconfounded treatment. While mflogit specifies the mean correctly, it is not robust to misspecification of the mean function. Therefore, wmflogit is the estimator of choice and is an appropriate doubly robust estimator for multivariate fractional outcome means. The only exception is when both the mean and PS models are misspecified, in which case all estimators are equally inconsistent in theory.

6.1. Population

First, we consider two covariates $(X_1,X_2^{\ast })$ where X₁ is continuously distributed, whereas X₂ is binary $$\begin{array}{l}{X}_2=\{\begin{array}{l}1,\text{ if }{X}_2^{\ast }+\frac{V}{3}>0\\ 0,\text{ otherwise }\end{array}\end{array}$$ where $V\sim N(0,1)$ and $(X_1,X_2^{\ast })$ is distributed bivariate normal with zero mean and variance-covariance matrix $$\begin{array}{l}{\mathrm{\Omega }}_X=\left(\begin{array}{cc}1& 0.4\\ 0.4& 1\end{array}\right)\end{array}$$

To generate the vector of treatment assignments, 𝐖_i, we use the following rule where the true propensity score is a multinomial logit with an index that includes square and interaction terms. $$\begin{array}{l}{\rho }_g(\mathbf{X})=\frac{\text{exp}(\mathbf{Z}{\mathit{\delta }}_g)}{\sum _{h=1}^3\text{exp}(\mathbf{Z}{\mathit{\delta }}_h)}\text{ for each }g=1,2,3\end{array}$$ where $\mathbf{Z}=(1,X_1,X_2,X_1{X}_2,X_1^2)$ and ${\mathit{\delta }}_g=({\delta }_{0g},{\delta }_{1g},{\delta }_{2g},{\delta }_{3g},{\delta }_{4g}{)}^{\prime }$. Using these probabilities and a uniform random variable, $U=\text{Uniform}(0,1)$, $$\begin{array}{l}{\mathbf{W}}_i=\{\begin{array}{l}(1,0,0),\text{ if }{U}_i\le {\rho }_1({\mathbf{X}}_i)\\ (0,1,0),\text{ if }{\rho }_1({\mathbf{X}}_i)<U_i\le {\rho }_1({\mathbf{X}}_i)+{\rho }_2({\mathbf{X}}_i)\\ (0,0,1),\text{ if }{\rho }_1({\mathbf{X}}_i)+{\rho }_2({\mathbf{X}}_i)<U_i\le 1.\end{array}\end{array}$$

The outcome of shares 𝐘(g) for each of the three treatment levels follows Dirichlet distribution denoted by Dir(𝜼_g) with parameters ${\mathit{\eta }}_g=({\eta }_{1g},{\eta }_{2g},{\eta }_{3g},{\eta }_{4g}{)}^{\prime }$ with ${\eta }_{jg}\equiv {\varphi }_g\cdot {\mu }_{jg}(\mathbf{X})$ and ${\eta }_{0g}\equiv \sum _{j=1}^4{\eta }_{jg}$ where we choose ϕ_g = 10. Given this parametrization, (21) $$\begin{array}{l}\mathbb{E}[Y_j(g)|\mathbf{X}]\equiv \frac{{\eta }_{jg}}{{\eta }_{0g}}={\mu }_{jg}(\mathbf{X})\end{array}$$(21) where $$\begin{array}{l}\begin{array}{c}{\mu }_{jg}(\mathbf{X})=\frac{\text{exp}(\mathbf{Z}{\mathit{\gamma }}_{jg})}{\sum _{h=1}^4\text{exp}(\mathbf{Z}{\mathit{\gamma }}_{jh})}\text{ for }j=1,2,3,4\text{ and }g=1,2,3\end{array}\end{array}$$ and ${\mathit{\gamma }}_{jg}=({\gamma }_{0jg},{\gamma }_{1jg},{\gamma }_{2jg},{\gamma }_{3jg},{\gamma }_{4jg})$. The parameter values for 𝜹_g and 𝜸_jg are given in Tables 2 and 3, respectively. Note that we normalize the outside option to 1 by setting the parameter values to 0. When estimating the wmflogit and mflogit estimators, misspecification of the mean and weights in Table 1 corresponds to mis-specifying the index, 𝐙𝜷, which includes squares and interaction terms to instead be linear in X₁ and X₂. Meanwhile, the multinomial functional form for the PO mean and propensity score is kept intact. Also note that, given this DGP, the linear means are always misspecified.

Table 2. Parameter values, δ, indexing the PS model.

g	δ0g	δ1g	δ2g	δ3g	δ4g
1	−0.2838	−0.2439	−0.1120	0.1077	−0.2119
2	0.1197	−0.1346	−0.2886	0.0820	−0.2173
3	0.0000	0.0000	0.0000	0.0000	0.0000

Table 3. Parameter values, γ, indexing the conditional mean model.

g	j	γ0jg	γ1jg	γ2jg	γ3jg	γ4jg
1	1	−0.8367	−0.5138	−0.7106	0.8645	0.7805
	2	0.5660	0.5274	−0.5115	0.1990	0.9790
	3	−0.0186	−0.8907	0.1897	0.6793	−0.9069
	4	0.0000	0.0000	0.0000	0.0000	0.0000
2	1	−0.7631	−0.6397	0.2590	0.3708	−0.8409
	2	−0.7480	−0.5753	0.5995	−0.7270	−0.8620
	3	−0.2406	0.2353	0.6972	0.5462	0.5900
	4	0.0000	0.0000	0.0000	0.0000	0.0000
3	1	−0.0750	−0.1683	−0.5333	0.8387	0.1446
	2	0.8224	−0.1144	−0.5736	0.0993	−0.2834
	3	−0.5627	0.3426	−0.7047	−0.6471	0.2987
	4	0.0000	0.0000	0.0000	0.0000	0.0000

6.2. Discussion

Tables 4–7 report bias and RMSE of lra, ipw, ipwlra, mflogit, and wmflogit for μ_jg, $j=1,2,3,4$ and $g=1,2,3$ for the four estimation scenarios listed in Table 1. The results from the simulation exercise line up with the theory. One can see that lra remains inconsistent for share means in all four cases. This bias does not diminish as sample size increases. For the cases where the weight is correctly specified (see Tables 4 and 5), all weighted estimators are consistent, whereas all unweighted estimators are inconsistent. For the case where the weight is misspecified but mean is correct (Table 6), estimators that depend crucially on the PS such as ipw and ipwlra are inconsistent, whereas wmflogit which is doubly robust remains consistent. Finally, case 4 corresponds to all elements of the proposed estimator being misspecified. Despite being completely wrong, we do see that wmflogit has the lowest RMSE relative to all other estimators for most treatment levels g and categories j even though all estimators are inconsistent in theory. Therefore, this article advocates for wmflogit which is doubly robust for multivariate fractional outcomes and exhibits favorable finite sample behavior.

Table 4. Bias and RMSE for μ_jg when mean and weights are both correct.

g	N=5,000
	j	lra		ipw		ipwlra		mflogit		wmflogit
		bias	rmse	bias	rmse	bias	rmse	bias	rmse	bias	rmse
1	1	0.0008	0.0014	−0.0003	0.0012	−0.0002	0.0012	−0.0003	0.0012	−0.0003	0.0012
	2	0.0038	0.0044	0.0001	0.0019	0.0006	0.0023	0.0001	0.0018	0.0001	0.0018
	3	−0.0031	0.0035	0.0000	0.0014	−0.0004	0.0017	0.0000	0.0012	0.0000	0.0012
	4	−0.0015	0.0020	0.0002	0.0012	−0.0001	0.0014	0.0002	0.0012	0.0002	0.0012
2	1	−0.0017	0.0019	−0.0001	0.0008	−0.0003	0.0009	0.0000	0.0008	0.0000	0.0008
	2	−0.0018	0.0021	0.0001	0.0011	0.0000	0.0012	0.0000	0.0010	0.0000	0.0010
	3	0.0041	0.0046	−0.0003	0.0019	0.0001	0.0023	−0.0002	0.0018	−0.0002	0.0018
	4	−0.0006	0.0014	0.0003	0.0012	0.0002	0.0013	0.0002	0.0012	0.0002	0.0012
3	1	0.0006	0.0013	0.0000	0.0012	0.0000	0.0012	0.0001	0.0012	0.0001	0.0012
	2	−0.0017	0.0021	0.0001	0.0013	0.0001	0.0012	0.0000	0.0012	0.0000	0.0012
	3	0.0007	0.0011	0.0000	0.0009	0.0000	0.0009	0.0000	0.0009	0.0000	0.0009
	4	0.0003	0.0011	−0.0001	0.0010	−0.0001	0.0010	−0.0001	0.0010	−0.0001	0.0010
N=8,000
1	1	0.0007	0.0010	−0.0003	0.0008	−0.0002	0.0008	−0.0003	0.0007	−0.0003	0.0007
	2	0.0042	0.0044	0.0005	0.0013	0.0007	0.0017	0.0004	0.0012	0.0004	0.0012
	3	−0.0032	0.0033	−0.0001	0.0008	−0.0003	0.0011	−0.0001	0.0008	−0.0001	0.0008
	4	−0.0017	0.0019	−0.0001	0.0007	−0.0002	0.0009	0.0000	0.0007	0.0000	0.0007
2	1	−0.0017	0.0018	−0.0001	0.0005	−0.0002	0.0006	−0.0001	0.0005	−0.0001	0.0005
	2	−0.0019	0.0020	0.0000	0.0006	−0.0001	0.0007	−0.0001	0.0006	−0.0001	0.0006
	3	0.0043	0.0045	0.0000	0.0011	0.0002	0.0014	0.0001	0.0011	0.0001	0.0011
	4	−0.0007	0.0011	0.0001	0.0008	0.0001	0.0008	0.0000	0.0008	0.0000	0.0007
3	1	0.0006	0.0010	0.0000	0.0008	0.0000	0.0008	0.0001	0.0008	0.0001	0.0008
	2	−0.0017	0.0019	0.0001	0.0008	0.0000	0.0008	0.0000	0.0008	0.0000	0.0008
	3	0.0008	0.0010	0.0000	0.0005	0.0000	0.0006	0.0000	0.0005	0.0000	0.0005
	4	0.0004	0.0007	−0.0001	0.0006	−0.0001	0.0006	−0.0001	0.0006	−0.0001	0.0006

^alra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

^bEmpirical distributions are generated with 1,000 Monte Carlo repetitions.

^cThe true population mean vector is: ${\mathit{\mu }}_1=(0.1245,0.5233,0.1648,0.1874),{\mathit{\mu }}_2=(0.0997,0.1305,0.5110,0.2588),{\mathit{\mu }}_3=(0.2550,0.3560,0.1342,0.2548)$.

Table 5. Bias and RMSE for μ_jg when mean is misspecified but weight is correct.

g	N=5,000
	j	lra		ipw		ipwlra		mflogit		wmflogit
		bias	rmse	bias	rmse	bias	rmse	bias	rmse	bias	rmse
1	1	−0.0042	0.0043	−0.0003	0.0012	−0.0003	0.0012	−0.0014	0.0018	−0.0002	0.0012
	2	−0.0142	0.0144	0.0001	0.0019	0.0001	0.0019	−0.0150	0.0152	0.0001	0.0019
	3	0.0097	0.0098	0.0000	0.0014	0.0000	0.0014	0.0089	0.0090	0.0000	0.0014
	4	0.0087	0.0088	0.0002	0.0012	0.0002	0.0012	0.0076	0.0077	0.0002	0.0012
2	1	0.0052	0.0053	−0.0001	0.0008	−0.0001	0.0008	0.0052	0.0053	−0.0001	0.0008
	2	0.0058	0.0059	0.0001	0.0011	0.0001	0.0011	0.0067	0.0068	0.0001	0.0011
	3	−0.0191	0.0193	−0.0003	0.0019	−0.0004	0.0020	−0.0192	0.0194	−0.0004	0.0019
	4	0.0081	0.0082	0.0003	0.0012	0.0004	0.0013	0.0073	0.0075	0.0003	0.0012
3	1	0.0111	0.0112	0.0000	0.0012	0.0000	0.0012	0.0103	0.0104	0.0000	0.0012
	2	−0.0150	0.0151	0.0001	0.0013	0.0001	0.0013	−0.0148	0.0149	0.0001	0.0013
	3	0.0075	0.0076	0.0000	0.0009	0.0000	0.0009	0.0077	0.0078	0.0000	0.0009
	4	−0.0036	0.0038	−0.0001	0.0010	−0.0001	0.0010	−0.0032	0.0034	−0.0001	0.0010
N=8,000
1	1	−0.0042	0.0043	−0.0003	0.0008	−0.0003	0.0007	−0.0014	0.0016	−0.0003	0.0007
	2	−0.0140	0.0140	0.0005	0.0013	0.0004	0.0013	−0.0148	0.0149	0.0004	0.0013
	3	0.0097	0.0097	−0.0001	0.0008	−0.0001	0.0009	0.0089	0.0089	−0.0001	0.0009
	4	0.0085	0.0086	−0.0001	0.0007	0.0000	0.0007	0.0074	0.0074	−0.0001	0.0008
2	1	0.0052	0.0053	−0.0001	0.0005	−0.0001	0.0005	0.0052	0.0052	−0.0001	0.0005
	2	0.0057	0.0057	0.0000	0.0006	0.0000	0.0007	0.0066	0.0067	0.0000	0.0007
	3	−0.0189	0.0189	0.0000	0.0011	−0.0001	0.0012	−0.0190	0.0190	0.0000	0.0012
	4	0.0079	0.0080	0.0001	0.0008	0.0002	0.0008	0.0072	0.0073	0.0001	0.0008
3	1	0.0111	0.0112	0.0000	0.0008	0.0000	0.0008	0.0103	0.0104	0.0000	0.0008
	2	−0.0151	0.0152	0.0001	0.0008	0.0001	0.0008	−0.0149	0.0149	0.0001	0.0008
	3	0.0076	0.0077	0.0000	0.0005	0.0000	0.0005	0.0078	0.0079	0.0000	0.0005
	4	−0.0036	0.0037	−0.0001	0.0006	−0.0001	0.0006	−0.0032	0.0033	−0.0001	0.0006

^alra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

^bEmpirical distributions are generated with 1,000 Monte Carlo repetitions.

^cThe true population mean vector is: ${\mathit{\mu }}_1=(0.1245,0.5233,0.1648,0.1874),{\mathit{\mu }}_2=(0.0997,0.1305,0.5110,0.2588),{\mathit{\mu }}_3=(0.2550,0.3560,0.1342,0.2548)$.

Table 6. Bias and RMSE for μ_jg when mean is correct but weight is misspecifed.

g	N=5,000
	j	lra		ipw		ipwlra		mflogit		wmflogit
		bias	rmse	bias	rmse	bias	rmse	bias	rmse	bias	rmse
1	1	0.0008	0.0014	−0.0026	0.0028	0.0006	0.0014	−0.0003	0.0012	−0.0003	0.0012
	2	0.0038	0.0044	−0.0149	0.0151	0.0033	0.0039	0.0001	0.0018	0.0001	0.0018
	3	−0.0031	0.0035	0.0092	0.0094	−0.0026	0.0030	0.0000	0.0012	0.0000	0.0012
	4	−0.0015	0.0020	0.0082	0.0083	−0.0013	0.0019	0.0002	0.0012	0.0002	0.0012
2	1	−0.0017	0.0019	0.0052	0.0053	−0.0015	0.0017	0.0000	0.0008	0.0000	0.0008
	2	−0.0018	0.0021	0.0058	0.0059	−0.0014	0.0018	0.0000	0.0010	0.0000	0.0010
	3	0.0041	0.0046	−0.0189	0.0190	0.0033	0.0039	−0.0002	0.0018	−0.0002	0.0018
	4	−0.0006	0.0014	0.0079	0.0080	−0.0004	0.0013	0.0002	0.0012	0.0002	0.0012
3	1	0.0006	0.0013	0.0117	0.0118	0.0007	0.0014	0.0001	0.0012	0.0001	0.0012
	2	−0.0017	0.0021	−0.0156	0.0157	−0.0014	0.0019	0.0000	0.0012	0.0000	0.0012
	3	0.0007	0.0011	0.0077	0.0078	0.0005	0.0010	0.0000	0.0009	0.0000	0.0009
	4	0.0003	0.0011	−0.0038	0.0040	0.0002	0.0010	−0.0001	0.0010	−0.0001	0.0010
N=8,000
1	1	0.0007	0.0010	−0.0026	0.0027	0.0006	0.0009	−0.0003	0.0007	−0.0003	0.0007
	2	0.0042	0.0044	−0.0147	0.0148	0.0036	0.0038	0.0004	0.0012	0.0004	0.0012
	3	−0.0032	0.0033	0.0092	0.0093	−0.0026	0.0028	−0.0001	0.0008	−0.0001	0.0008
	4	−0.0017	0.0019	0.0080	0.0081	−0.0015	0.0018	0.0000	0.0007	0.0000	0.0007
2	1	−0.0017	0.0018	0.0052	0.0052	−0.0015	0.0016	−0.0001	0.0005	−0.0001	0.0005
	2	−0.0019	0.0020	0.0057	0.0057	−0.0015	0.0016	−0.0001	0.0006	−0.0001	0.0006
	3	0.0043	0.0045	−0.0186	0.0187	0.0035	0.0038	0.0001	0.0011	0.0001	0.0011
	4	−0.0007	0.0011	0.0078	0.0078	−0.0006	0.0010	0.0000	0.0008	0.0000	0.0007
3	1	0.0006	0.0010	0.0118	0.0118	0.0007	0.0010	0.0001	0.0008	0.0001	0.0008
	2	−0.0017	0.0019	−0.0158	0.0158	−0.0015	0.0017	0.0000	0.0008	0.0000	0.0008
	3	0.0008	0.0010	0.0078	0.0078	0.0005	0.0008	0.0000	0.0005	0.0000	0.0005
	4	0.0004	0.0007	−0.0038	0.0039	0.0003	0.0007	−0.0001	0.0006	−0.0001	0.0006

^alra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

^bEmpirical distributions are generated with 1,000 Monte Carlo repetitions.

^cThe true population mean vector is: ${\mathit{\mu }}_1=(0.1245,0.5233,0.1648,0.1874),{\mathit{\mu }}_2=(0.0997,0.1305,0.5110,0.2588),{\mathit{\mu }}_3=(0.2550,0.3560,0.1342,0.2548)$.

Table 7. Bias and RMSE for μ_jg when mean and weight are both misspecified.

g	N=5,000
	j	lra		ipw		ipwlra		mflogit		wmflogit
		bias	rmse	bias	rmse	bias	rmse	bias	rmse	bias	rmse
1	1	−0.0042	0.0043	−0.0026	0.0028	−0.0037	0.0039	−0.0014	0.0018	−0.0014	0.0019
	2	−0.0142	0.0144	−0.0149	0.0151	−0.0128	0.0130	−0.0150	0.0152	−0.0134	0.0136
	3	0.0097	0.0098	0.0092	0.0094	0.0087	0.0088	0.0089	0.0090	0.0081	0.0083
	4	0.0087	0.0088	0.0082	0.0083	0.0078	0.0079	0.0076	0.0077	0.0067	0.0069
2	1	0.0052	0.0053	0.0052	0.0053	0.0050	0.0051	0.0052	0.0053	0.0049	0.0050
	2	0.0058	0.0059	0.0058	0.0059	0.0054	0.0056	0.0067	0.0068	0.0065	0.0067
	3	−0.0191	0.0193	−0.0189	0.0190	−0.0182	0.0184	−0.0192	0.0194	−0.0184	0.0185
	4	0.0081	0.0082	0.0079	0.0080	0.0077	0.0078	0.0073	0.0075	0.0069	0.0070
3	1	0.0111	0.0112	0.0117	0.0118	0.0125	0.0126	0.0103	0.0104	0.0121	0.0122
	2	−0.0150	0.0151	−0.0156	0.0157	−0.0160	0.0161	−0.0148	0.0149	−0.0159	0.0159
	3	0.0075	0.0076	0.0077	0.0078	0.0075	0.0076	0.0077	0.0078	0.0076	0.0077
	4	−0.0036	0.0038	−0.0038	0.0040	−0.0040	0.0041	−0.0032	0.0034	−0.0038	0.0039
N=8,000
1	1	−0.0042	0.0043	−0.0026	0.0027	−0.0037	0.0038	−0.0014	0.0016	−0.0014	0.0016
	2	−0.0140	0.0140	−0.0147	0.0148	−0.0126	0.0126	−0.0148	0.0149	−0.0132	0.0133
	3	0.0097	0.0097	0.0092	0.0093	0.0086	0.0087	0.0089	0.0089	0.0081	0.0081
	4	0.0085	0.0086	0.0080	0.0081	0.0076	0.0077	0.0074	0.0074	0.0066	0.0066
2	1	0.0052	0.0053	0.0052	0.0052	0.0050	0.0051	0.0052	0.0052	0.0049	0.0050
	2	0.0057	0.0057	0.0057	0.0057	0.0053	0.0054	0.0066	0.0067	0.0064	0.0065
	3	−0.0189	0.0189	−0.0186	0.0187	−0.0179	0.0180	−0.0190	0.0190	−0.0181	0.0182
	4	0.0079	0.0080	0.0078	0.0078	0.0076	0.0076	0.0072	0.0073	0.0067	0.0068
3	1	0.0111	0.0112	0.0118	0.0118	0.0126	0.0126	0.0103	0.0104	0.0121	0.0122
	2	−0.0151	0.0152	−0.0158	0.0158	−0.0162	0.0162	−0.0149	0.0149	−0.0160	0.0160
	3	0.0076	0.0077	0.0078	0.0078	0.0076	0.0076	0.0078	0.0079	0.0076	0.0077
	4	−0.0036	0.0037	−0.0038	0.0039	−0.0040	0.0041	−0.0032	0.0033	−0.0038	0.0038

^alra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

^bEmpirical distributions are generated with 1,000 Monte Carlo repetitions.

^cThe true population mean vector is: ${\mathit{\mu }}_1=(0.1245,0.5233,0.1648,0.1874),{\mathit{\mu }}_2=(0.0997,0.1305,0.5110,0.2588),{\mathit{\mu }}_3=(0.2550,0.3560,0.1342,0.2548)$.

7. Women’s time-use with Progresa

In this section, we use data from the Mexican conditional cash transfer program, Progresa, to construct time-use shares for women aged 15-65 years corresponding to three broad types of activities, namely, domestic work, farm work, and market work. Our objective here is to estimate the mean proportion of time spent by women using wmflogit in each of these activities depending upon their participation status in the program.

Progresa is a well-studied conditional cash transfer program which began in Mexico in 1997. The objective of the program was to encourage investment in human capital and reduce the transmission of poverty by offering cash assistance to women in poor Mexican households conditional on children’s school enrollment, school attendance, and visits to health clinics.

Eligibility for the program was determined in two phases where the first phase of the study used a deprivation index to determine which localities or communities in Mexico were highly deprived. A second phase of the program collected data on socioeconomic characteristics of households identified in the first phase using the Survey of Household Socio-Economic Characteristics (ENCASEH survey). Each locality selected in this phase was then randomly allocated to the intervention or control group. A subsequent poverty index was used to determine eligibility at the household level and those that were eligible received program benefits in the intervention localities.

We use the Evaluation Survey of Progresa (ENCEL surveys) which is a panel dataset with information on all rounds of data collection before and after Progresa was implemented. This panel is used to identify the baseline survey (round one) along with the time-use survey (round seven) to construct the final dataset used in this empirical study. Even though the original sample contains 24,077 households (142,697 observations), we restrict our analysis to households that have complete information on time-use shares and covariates for women between ages 15-65 years. This helps to ensure comparability of estimates across the different regression specifications.

We use this sample to study the effect of being a beneficiary or non beneficiary household on the time-use of women. To construct the dummy variable indicating the participation status of a given household, we interact the treatment-control locality dummy, contba, along with the dummy for being poor, pobre. This is important since only eligible households in the intervention localities were allowed to be Progresa participants.

The time-use shares for these women are then calculated by classifying 16 different types of activities (see Table 8) recorded in the time-use module of Progresa into three broad categories, namely, domestic work, farm work, and market work.⁸ Shares are then calculated by dividing time spent on each category by the total time spent in all categories.

Table 8. Types of activity in Progresa’s time-use module.

	Type of activity
1	Carrying water
2	Caring for or raising household animals
3	Cooking or preparing breakfast
4	Grocery shopping
5	Performing a free service
6	Washing and ironing clothes
7	Cleaning the house
8	Dropping off or picking up a member
9	Working on own business
10	Other activities
11	Making repairs to home
12	Knitting
13	Working the land on the plot
14	Working for a wage
15	Commuting to work
16	Performing a free service

Note that while the activity “Performing a free service” appears twice in the above list, one pertains to delivering a free service at home and the other at work.

The share equation in each type of regression controls for both individual level characteristics such as age, education, relationship with the head, and household level factors like household size, assets (such as owning a television, radio, refrigerator and other consumer durables), availability of electricity, poverty index, total agricultural land, number of dependents, head’s gender and education. The PS is estimated as a logit which controls for the marginalization index⁹ and characteristics that were used to determine households’ poverty status in order to be eligible to receive program benefits. The latter includes running water, electricity, ownership of durable goods, animals and land, number of dependents, and presence of disabled individuals (Parker and Skoufias, 2000; Schultz, 2004).

We then use lra, ipw, ipwlra, mflogit, and wmflogit to estimate PO means for beneficiary and non beneficiary households (see Table 9). The proposed estimator is implemented by using Stata’s fmlogit command developed by Buis (2008). These PO means are then used to estimate the effect of Progresa on time spent by women in each category of work. Table 10 reports the average treatment effect estimates for Progressa beneficiaries using the proposed wmflogit and the other above mentioned alternatives. We find that the difference in time spent on domestic work between beneficiary and non beneficiary households is marginally insignificant when using wmflogit and ipwlra but significant if one uses ipw and mflogit. Similary, the difference in time spent on market work is highly insignificant when using wmflogit and ipwlra but significant if one believes lra, ipw, and mflogit. However, for the case of farm activity, the difference is highly significant with wmflogit and insignificant if one believes the other estimators. The fact that wmflogit is doubly robust should engender more confidence in the results obtained from it compared to the other estimators considered.

Table 9. Estimate of potential time-use share means (PO means) by intervention group.

Progresa membership	Time-use shares	lra	ipw	ipwlra	mflogit	wmflogit
Intended non beneficiaries	domestic	0.7825	0.7853	0.7830	0.7822	0.7777
		(0.0037)	(0.0039)	(0.0038)	(0.0035)	(0.0043)
	farm	0.0202	0.0213	0.0207	0.0206	0.0187
		(0.0011)	(0.0012)	(0.0011)	(0.0011)	(0.0012)
	market	0.1972	0.1934	0.1963	0.1972	0.2037
		(0.0036)	(0.0038)	(0.0037)	(0.0034)	(0.0042)
Intended beneficiaries	domestic	0.7691	0.7458	0.7712	0.7680	0.8052
		(0.0062)	(0.0083)	(0.0069)	(0.0056)	(0.0127)
	farm	0.0188	0.0207	0.0197	0.0163	0.0131
		(0.0020)	(0.0030)	(0.0027)	(0.0021)	(0.0025)
	market	0.2121	0.2335	0.2091	0.2157	0.1817
		(0.0060)	(0.0080)	(0.0064)	(0.0053)	(0.0122)
Observations		15,777	15,777	15,777	15,777	15,777

^alra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

^bStandard error in parentheses.

^cThe share equations control for the woman’s age, education, relationship to the household head (wife or child), durable assets, availability of electricity, agricultural vehicles (like tractors, threshing machines etc), total agricultural land (in hectares), total number of dependents, household size, household head’s age, sex, and poverty index of the household.

^dPS is estimated as a logit with controls such as marginalization index, availability of running water, total number of dependents, ownership of consumer durables, animals, agricultural land, and presence of disabled individuals.

Table 10. Average treatment effect estimates.

Time-use shares	lra	ipw	ipwlra	mflogit	wmflogit
domestic	−0.0135	−0.0395	−0.0118	−0.0139	0.0284
	(0.0070)	(0.0091)	(0.0077)	(0.0067)	(0.0147)
	[0.0540]	[0.0000]	[0.1230]	[0.0390]	[0.0530]
farm	−0.0014	−0.0006	−0.0010	−0.0045	−0.0057
	(0.0023)	(0.0032)	(0.0029)	(0.0023)	(0.0012)
	[0.5420]	[0.8550]	[0.7330]	[0.0540]	[0.0000]
market	0.0149	0.0401	0.0128	0.0184	−0.0227
	(0.0067)	(0.0088)	(0.0072)	(0.0065)	(0.0143)
	[0.0270]	[0.0000]	[0.0730]	[0.0050]	[0.1130]

^a.lra refers to linear regression adjustment, ipw refers to inverse probability weighted sample averages for the different treatment levels, ipwlra refers to weighted linear regression adjustment, mflogit refers to multinomial fractional logit and finally wmflogit is the proposed weighted multinomial fractional logit estimator.

^b.Standard error in parentheses. The mflogit and wmflogit standard error values have been bootstrapped using 1000 replications.

^c.p-values are given in square brackets.

^d.The share equations control for the woman’s age, education, relationship to the household head (wife or child), durable assets, availability of electricity, agricultural vehicles (like tractors, threshing machines etc), total agricultural land (in hectares), total number of dependents, household size, household head’s age, sex, and poverty index of the household.

^e.PS is estimated as a logit with controls such as marginalization index, availability of running water, total number of dependents, ownership of consumer durables, animals, agricultural land, and presence of disabled individuals.

8. Conclusion

This article proposes a doubly robust estimator for estimating the vector of potential outcome means for outcomes that are expressed as a vector of mutually exclusive and exhaustive binary outcomes, shares, or proportions. Our framework allows for multivalued treatments, which means that all pairwise treatments effects are also estimated in a doubly robust manner. The estimator is a propensity score weighted multivariate fractional logit, termed wmflogit, which uses features of the linear exponential family of distributions to guarantee robustness to misspecification of either the propensity score model or the multivariate fractional logit mean as long as both are not misspecified at the same time. A practical advantage of this estimator is that it only requires one of the features (either the conditional mean or the propensity score) to be correctly specified thereby offering practitioners some hope of obtaining valid treatment effect estimates even in the face of uncertainty about the accuracy of model specification.

Our simulations illustrate the double robustness property quite well. We generate a 4×1 vector of shares from a Dirichlet distribution for three treatment levels. On comparing wmflogit with linear, weighted linear, inverse probability weighted, and unweighted multivariate fractional logit (mflogit) regressions, we find that the RMSE associated with wmflogit declines as sample size increases in all cases including the one where mean and weights are both incorrect. For mflogit, this happens only when the mean model is correctly specified. Linear regression is never consistent and that is easily seen in the simulations where the finite sample bias persists even as the sample size increases. Weighted linear regression adjustment (ipwlra) has an RMSE comparable to wmflogit only when the weight is correct. Since the mean will rarely be well represented by a linear function for multivariate fractional outcomes, the consistency argument for ipwlra rests entirely on correct weight. A similar argument holds for ipw. The only procedure that has a real chance at being doubly robust is wmflogit and this is easily observed in the simulations.

An application to time-use shares of women from the Mexican conditional cash transfer program, Progresa, also shows that the significance of ATE estimates is different when using wmflogit vs. when using other linear and unweighted alternatives. Given that wmflogit is robust to inconsistency in two directions, it should invoke more confidence in the results obtained from wmflogit compared to other alternatives when explaining multivariate fractional outcomes.

Acknowledgments

We would like to thank Dalia Ghanem and Karen del Mar Ortiz Becerra for providing us the Progresa dataset.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1. A version, known as AIPW, was proposed and studied in a string of papers by Robins et al. (1994), Robins et al. (2000), and Scharfstein et al. (1999). More recently, improved versions of the AIPW estimator have been proposed in both the missing data and causal inference literatures.

2. Including pre-treatment outcomes in the covariate vector makes unconfoundedness even more plausible as an identifying assumption. As in the case of Hirano and Imbens (2001), unconfoundedness is also easy to justify when we have a rich set of controls in the data that make it approximately true.

3. This QLL can be obtained by applying the general result on weighting given in Lemma 3.2 of Sloczyński and Wooldridge (2018) to the multinomial QLL function.

4. As mentioned in Imbens (2000), one may estimate the propensity scores using discrete response models if there is no natural ordering among the treatment levels or ordered response models if there is a natural ordering to the alternatives.

5. See section A for a derivation of the first-order conditions for the wmflogit estimator.

6. Similar equations for the ATE and the ATT can also be found in Hirano et al. (2003) for the binary case.

7. These sample sizes ensure that we have enough variation in the different share categories across the three treatment levels.

8. Two school related activities, namely, min_esc and min_tar are not used in construction of any of the three activities.

9. Progresa used marginalization index to geographically target highly marginalized localities in the seven states of Mexico. These localities were then randomized to either be in the intervention or control group.

Appendix A

first order conditions for weighted multivariate fractional logit

The first-order conditions for the weighted multinomial logistic regression can be derived by considering the following quasi maximum likelihood problem $$\begin{array}{l}L={\mathrm{\Pi }}_{i=1}^N{\mathrm{\Pi }}_{j=1}^J\{\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot [m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g){]}^{Y_{ij}}\}\end{array}$$

Then, the multinomial log likelihood function is given as (A.1) $$\begin{array}{l}\begin{array}{ll}\text{ln(L)}& =\sum _{i=1}^N{l}_i({\hat{\mathit{\gamma }}}_g,\hat{\mathit{\delta }})\text{ where }{l}_i({\hat{\mathit{\gamma }}}_g,\hat{\mathit{\delta }})=\sum _{j=1}^J\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot Y_{ij}\cdot \text{log}[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)]\end{array}\end{array}$$(A.1)

Now differentiating the above with respect to ${\hat{\mathit{\gamma }}}_g$ gives us the following vector of first-order conditions $$\begin{array}{ll}\sum _{i=1}^N{\mathbf{s}}_{i1}({\hat{\mathit{\gamma }}}_g,\hat{\mathit{\delta }})& ={\mathbf{0}}_{(J-1)(K+1)\times 1}\end{array}$$ where $s_{i1}({\mathit{\gamma }}_g,\widehat{\mathit{\delta }})\equiv \nabla {\mathit{\gamma }}_g{l}_i({\mathit{\gamma }}_g,\widehat{\mathit{\delta }})=\left(\begin{array}{c}{\frac{\partial l_i}{\partial {\mathit{\gamma }}_{1g}}}^{\prime }\\ ⋮\\ {\frac{\partial l_i}{\partial {\mathit{\gamma }}_{jg}}}^{\prime }\\ ⋮\\ {\frac{\partial l_i}{\partial {\mathit{\gamma }}_{(J-1g)}}}^{\prime }\end{array}\right)=\left(\begin{array}{c}{s}_{i11}\\ ⋮\\ s_{i1j}\\ ⋮\\ s_{i1(J-1)}\end{array}\right)$ and for each $j=1,\dots ,J-1$ (A.2) $$\begin{array}{l}{\mathbf{s}}_{i1j}=\left(\begin{array}{c}\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot [Y_{ij}-m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]\\ \frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot {\mathbf{X}}_i^{\prime }\cdot [Y_{ij}-m_{jg}({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]\end{array}\right)={\mathbf{0}}_{(K+1)\times 1}.\end{array}$$(A.2)

If we stack first order conditions corresponding to all the intercepts and slopes separately, one may write them compactly as (A.3) $$\begin{array}{ll}\sum _{i=1}^N\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot [{\mathbf{Y}}_i-{\mathbf{m}}_g({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]& ={\mathbf{0}}_{(J-1)\times 1}\end{array}$$(A.3) (A.4) $$\begin{array}{ll}\sum _{i=1}^N\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot {\mathbf{X}}_i^{\prime }\cdot [{\mathbf{Y}}_i-{\mathbf{m}}_g({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]& ={\mathbf{0}}_{K(J-1)\times 1}\end{array}$$(A.4)

Similarly, the first-order conditions corresponding to ${\hat{\mathit{\gamma }}}_g$ which solves (12(12) $$\begin{array}{l}\sum _{i=1}^N\sum _{j=1}^J\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}{p}_h({\mathbf{X}}_i,\hat{\mathit{\delta }})\cdot Y_{ij}\cdot \text{log}[m_{jg}({\mathbf{X}}_i,{\mathit{\gamma }}_g)].\end{array}$$(12) ) can be stacked to obtain (A.5) $$\begin{array}{ll}\sum _{i=1}^N\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot p_h({\mathbf{X}}_i,\hat{\mathit{\delta }})\cdot [{\mathbf{Y}}_i-{\mathbf{m}}_g({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]& ={\mathbf{0}}_{(J-1)\times 1}\end{array}$$(A.5) (A.6) $$\begin{array}{ll}\sum _{i=1}^N\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot p_h({\mathbf{X}}_i,\hat{\mathit{\delta }})\cdot {\mathbf{X}}_i^{\prime }\otimes [{\mathbf{Y}}_i-{\mathbf{m}}_g({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]& ={\mathbf{0}}_{K(J-1)\times 1}\end{array}$$(A.6)

Appendix B

Proofs

Proof of Proposition 1.

Case 1: When the propensity score model is correctly specified: In this case, the population analogue of the sample first-order conditions corresponding to the intercept in (A.3(A.3) $$\begin{array}{ll}\sum _{i=1}^N\frac{W_{ig}}{p_g({\mathbf{X}}_i,\hat{\mathit{\delta }})}\cdot [{\mathbf{Y}}_i-{\mathbf{m}}_g({\mathbf{X}}_i,{\hat{\mathit{\gamma }}}_g)]& ={\mathbf{0}}_{(J-1)\times 1}\end{array}$$(A.3) ) is given by (B.1) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot (\mathbf{Y}-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))]={\mathbf{0}}_{J\times 1}\end{array}$$(B.1) where ${\mathit{\gamma }}_g^{\ast }\ne {\mathit{\gamma }}_g^0$. The above can be rearranged to give us (B.2) $$\begin{array}{l}\begin{array}{c}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y})=\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))\end{array}\end{array}$$(B.2)

Using (2.3(3) $$\begin{array}{l}\mathbf{Y}=W_1\cdot \mathbf{Y}(1)+\dots +W_G\cdot \mathbf{Y}(G).\end{array}$$(3) ) we know that for each g, $$\begin{array}{l}\begin{array}{ll}{W}_g\cdot \mathbf{Y}& =W_g\cdot (W_0\mathbf{Y}(0)+W_1\mathbf{Y}(1)+\dots +W_G\mathbf{Y}(G))\\ & =W_g\cdot \mathbf{Y}(g)\end{array}\end{array}$$

Rewriting the left hand side of Eq. (B.2(B.2) $$\begin{array}{l}\begin{array}{c}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y})=\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))\end{array}\end{array}$$(B.2) ), $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y}]=\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y}(g)]\end{array}$$

Using law of iterated expectations, we can rewrite $$\begin{array}{l}\begin{array}{ll}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y}(g)]& =\mathbb{E}[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y}(g)|\mathbf{X},\mathbf{W})]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbb{E}(\mathbf{Y}(g)|\mathbf{X},\mathbf{W})]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbb{E}(\mathbf{Y}(g)|\mathbf{X})]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathit{\mu }}_g(\mathbf{X})]\end{array}\end{array}$$

Using iterated expectations again, $$\begin{array}{l}\begin{array}{ll}& =\mathbb{E}[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathit{\mu }}_g(\mathbf{X})|\mathbf{X})]\\ & =\mathbb{E}[\frac{{\mathit{\mu }}_g(\mathbf{X})}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbb{P}(W_g=1|\mathbf{X})]\\ & =\mathbb{E}[\mathbf{Y}(g)]\end{array}\end{array}$$

Similarly, rewrite the right hand side of Eq. (B.2(B.2) $$\begin{array}{l}\begin{array}{c}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{Y})=\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))\end{array}\end{array}$$(B.2) ), $$\begin{array}{l}\begin{array}{ll}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]& =\mathbb{E}[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|\mathbf{X})]\\ & =\mathbb{E}[\frac{{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbb{P}(W_g=1|\mathbf{X})]\\ & =\mathbb{E}[{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\end{array}\end{array}$$

Therefore, combining the two results we get (B.3) $$\begin{array}{l}\mathbb{E}[\mathbf{Y}(g)]=\mathbb{E}[{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\end{array}$$(B.3)

Part 2: For 𝝁_g|h, the population first-order conditions for the intercept are given by (B.4) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot (\mathbf{Y}-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))]={\mathbf{0}}_{J\times 1}\end{array}$$(B.4) which can again be re-arranged to give us $$\begin{array}{l}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{Y})=\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))\end{array}$$

Consider the left hand side, $$\begin{array}{l}\begin{array}{ll}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{Y})& =\mathbb{E}\left[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{Y}(g)|\mathbf{X},\mathbf{W})\right]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbb{E}(\mathbf{Y}(g)|\mathbf{X})]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathit{\mu }}_g(\mathbf{X})]\\ & =\mathbb{E}\left[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathit{\mu }}_g(\mathbf{X})|\mathbf{X})\right]\\ & =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_0)}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot {\mathit{\mu }}_g(\mathbf{X})\cdot \mathbb{P}(W_g=1|\mathbf{X})]\\ & =\mathbb{E}[p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathit{\mu }}_g(\mathbf{X})]\end{array}\end{array}$$

Now, (B.5) $$\begin{array}{l}\begin{array}{ll}\mathbb{E}[p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathit{\mu }}_g(\mathbf{X})]& =\mathbb{E}[\mathbb{P}(W_h=1|\mathbf{X})\cdot \mathbb{E}(\mathbf{Y}(g)|\mathbf{X})]\\ & =\mathbb{E}[W_h\cdot \mathbf{Y}(g)]\\ & =\mathbb{E}[\mathbf{Y}(g)|W_h=1]\cdot \mathbb{P}(W_h=1)\end{array}\end{array}$$(B.5)

Now let’s consider the right hand side, (B.6) $$\begin{array}{l}\begin{array}{ll}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))& =\mathbb{E}\left[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|\mathbf{Y}(g),\mathbf{X})\right]\\ & =\mathbb{E}[\frac{1}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })\cdot \mathbb{P}(W_g=1|\mathbf{Y}(g),\mathbf{X})]\\ & =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_0)}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })\cdot \mathbb{P}(W_g=1|\mathbf{X})]\\ & =\mathbb{E}[p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\\ & =\mathbb{E}[\mathbb{E}(W_h\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|\mathbf{X})]\\ & =\mathbb{E}[W_h\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\\ & =\mathbb{P}(W_h=1)\cdot \mathbb{E}[\mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|W_h=1]\end{array}\end{array}$$(B.6)

Now, comparing Eqs. (B.5(B.5) $$\begin{array}{l}\begin{array}{ll}\mathbb{E}[p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathit{\mu }}_g(\mathbf{X})]& =\mathbb{E}[\mathbb{P}(W_h=1|\mathbf{X})\cdot \mathbb{E}(\mathbf{Y}(g)|\mathbf{X})]\\ & =\mathbb{E}[W_h\cdot \mathbf{Y}(g)]\\ & =\mathbb{E}[\mathbf{Y}(g)|W_h=1]\cdot \mathbb{P}(W_h=1)\end{array}\end{array}$$(B.5) ) and (B.6(B.6) $$\begin{array}{l}\begin{array}{ll}\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast }))& =\mathbb{E}\left[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|\mathbf{Y}(g),\mathbf{X})\right]\\ & =\mathbb{E}[\frac{1}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })\cdot \mathbb{P}(W_g=1|\mathbf{Y}(g),\mathbf{X})]\\ & =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_0)}{p_g(\mathbf{X},{\mathit{\delta }}_0)}\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })\cdot \mathbb{P}(W_g=1|\mathbf{X})]\\ & =\mathbb{E}[p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\\ & =\mathbb{E}[\mathbb{E}(W_h\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|\mathbf{X})]\\ & =\mathbb{E}[W_h\cdot \mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })]\\ & =\mathbb{P}(W_h=1)\cdot \mathbb{E}[\mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|W_h=1]\end{array}\end{array}$$(B.6) ), $$\begin{array}{l}\mathbb{E}[\mathbf{m}(\mathbf{X},{\mathit{\gamma }}_g^{\ast })|W_h=1]={\mathit{\mu }}_{g|h}\end{array}$$

Case 2: When the mean model is correctly specified

In this case, the population first order conditions are given by (B.7) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot (\mathbf{Y}-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]={\mathbf{0}}_{J\times 1}\text{ for each }g=1,\dots ,G\end{array}$$(B.7) where ${\mathit{\delta }}_g^{\ast }$ may or may not be equal to the true propensity score parameter, 𝜹₀. Rewriting the left hand side of eq (B.7(B.7) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot (\mathbf{Y}-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]={\mathbf{0}}_{J\times 1}\text{ for each }g=1,\dots ,G\end{array}$$(B.7) ) and using the fact that for each g, $$\begin{array}{l}\begin{array}{ll}{W}_g\cdot \mathbf{Y}& =W_g\cdot (W_0\mathbf{Y}(0)+W_1\mathbf{Y}(1)+\dots +W_G\mathbf{Y}(G))\\ & =W_g\cdot \mathbf{Y}(g)\end{array}\end{array}$$ we get (B.8) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]=\mathbf{0}\end{array}$$(B.8)

Using law of iterated expectations, $$\begin{array}{l}\begin{array}{ll}\mathbb{E}[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))|\mathbf{X},\mathbf{W})]& =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{\mathbb{E}(\mathbf{Y}(g)|\mathbf{X},\mathbf{W})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{\mathbb{E}(\mathbf{Y}(g)|\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\\ & =\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\end{array}\end{array}$$

Using iterated expectations again, (B.9) $$\begin{array}{ll}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]& =\mathbb{E}[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}|\mathbf{X})]\\ & =\mathbb{E}[\frac{\{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \mathbb{P}(W_g=1|\mathbf{X})]\\ & =\mathbb{E}[\frac{{\rho }_g(\mathbf{X})}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\end{array}$$(B.9)

Therefore, combining Eq. (B.9(B.5) $$\begin{array}{l}\begin{array}{ll}\mathbb{E}[p_h(\mathbf{X},{\mathit{\delta }}_0)\cdot {\mathit{\mu }}_g(\mathbf{X})]& =\mathbb{E}[\mathbb{P}(W_h=1|\mathbf{X})\cdot \mathbb{E}(\mathbf{Y}(g)|\mathbf{X})]\\ & =\mathbb{E}[W_h\cdot \mathbf{Y}(g)]\\ & =\mathbb{E}[\mathbf{Y}(g)|W_h=1]\cdot \mathbb{P}(W_h=1)\end{array}\end{array}$$(B.5) ) with (B.8(B.8) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]=\mathbf{0}\end{array}$$(B.8) ), we get (B.10) $$\begin{array}{l}\mathbb{E}[\frac{p_g(\mathbf{X},{\mathit{\gamma }}_g^0)}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]=\mathbf{0}\end{array}$$(B.10)

Since $\frac{{\rho }_g(\mathbf{X})}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}>0$ the above will be true only when ${\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)={\mathit{\mu }}_g(\mathbf{X})$ or in order words, when the mflogit mean is correctly specified.

Part 2: For 𝝁_g|h, consider the population first-order conditions given by (B.11) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })\cdot (\mathbf{Y}-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]={\mathbf{0}}_{J\times 1}\end{array}$$(B.11)

Again, because of $W_g\cdot \mathbf{Y}=W_g\cdot \mathbf{Y}(g)$, (B.12) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]\\ =\mathbb{E}\left[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))|\mathbf{X},\mathbf{W})\right]\\ =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot W_g\cdot \{\mathbb{E}(\mathbf{Y}(g)|\mathbf{X},\mathbf{W})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\\ =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot W_g\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\\ =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot {\rho }_g(\mathbf{X})\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\end{array}$$(B.12) where the second equality uses iterated expectations, fourth equality uses unconfoundedness, and the last equality uses a second application of iteration expectations. Given that $\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot {\rho }_g(\mathbf{X})>0$, Eqs. (B.11(B.11) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })\cdot (\mathbf{Y}-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]={\mathbf{0}}_{J\times 1}\end{array}$$(B.11) ) and (B.12(B.12) $$\begin{array}{l}\mathbb{E}[\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))]\\ =\mathbb{E}\left[\mathbb{E}(\frac{W_g}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })\cdot (\mathbf{Y}(g)-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0))|\mathbf{X},\mathbf{W})\right]\\ =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot W_g\cdot \{\mathbb{E}(\mathbf{Y}(g)|\mathbf{X},\mathbf{W})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\\ =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot W_g\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\\ =\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot {\rho }_g(\mathbf{X})\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]\end{array}$$(B.12) ) implies that (B.13) $$\begin{array}{l}\mathbb{E}[\frac{p_h(\mathbf{X},{\mathit{\delta }}_h^{\ast })}{p_g(\mathbf{X},{\mathit{\delta }}_g^{\ast })}\cdot {\rho }_g(\mathbf{X})\cdot \{{\mathit{\mu }}_g(\mathbf{X})-{\mathbf{m}}_g(\mathbf{X},{\mathit{\gamma }}_g^0)\}]=\mathbf{0}\end{array}$$(B.13) only when the multinomial logit mean is correctly specified. ▪