4 Weighted log composite likelihood: a unified approach

J.N.K. Rao, F. Verret and M.A. Hidiroglou

Previous | Next

In this section we propose a unified approach applicable to both linear and generalized linear multi-level models. This approach is based on the concept of composite likelihood which has become popular in the non-survey literature to handle clustered or spatial data (see e.g., Lindsay 1988, Lele and Taper 2002 and Varin, Reid and Firth 2011). A pairwise marginal composite likelihood is obtained by multiplying the likelihood contributions from all the distinct pairs within clusters. Note that the composite likelihood is obtained by pretending the sub-models are independent. When the super-population model holds for the sample, then we can obtain parameter estimators by maximizing the pairwise composite likelihood. Here we extend this approach to handle informative designs by obtaining weighted estimating equations that require only the marginal weights $w_i$ and $w_{j|i}$ and the pairwise weights $w_{jk|i},$ as in Section 3.

The census log pairwise composite likelihood is given by

$$l_C(\theta )={\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j<k=1}^{M_i}\log f\left(y_{ij},y_{ik}|\theta \right),}}\left(4.1\right)$$

where $f\left(y_{ij},y_{ik}|\theta \right)$ is the marginal joint density of $y_{ij}$ and $y_{ik}.$ We estimate (4.1) by the design-weighted log pairwise composite likelihood

$$l_{wC}(\theta )={\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}\log f\left(y_{ij},y_{ik}|\theta \right)}}\left(4.2\right)$$

which depends only on the first order level 1 and level 2 inclusion probabilities and the second order level 1 probabilities. We then solve the weighted composite score equations

$${\widehat{U}}_{wC}(\theta )=\partial l_{wC}(\theta )/\partial \theta =0,\left(4.3\right)$$

obtained from (4.2) to get a weighted composite likelihood estimator, ${\widehat{\theta }}_{wC},$ of $\theta$ . The proposed method is applicable to linear and generalized linear two-level models.

We note that ${\widehat{U}}_{wC}(\theta )$ , given by (4.3), is a vector of estimating functions with zero expectation with respect to the design and the model, i.e., $E_m{E}_p\left\{{\widehat{U}}_{wC}(\theta )\right\}=0.$ Using this result, it can be shown that the weighted composite likelihood (WCL) estimator ${\widehat{\theta }}_{wC}$ of $\theta$ is design-model consistent as the number of level 2 units in the sample, $n,$ increases, even when the within cluster sample sizes, $m_i,$ are small. Details of the proof are given in Yi, Rao and Li (2012). In the non-survey context, we have limited theoretical and empirical evidence that the composite likelihood approach leads to efficient estimators (e.g., Bellio and Varin 2005, Lindsay et al. 2011). Our simulation study (Section 5) indicates that the weighted composite likelihood approach performs well in terms of efficiency, even for small within-cluster sample sizes.

In the case of the nested error model (3.13), following Lele and Taper (2002) we can simplify the pairwise composite likelihood approach by replacing the bivariate density function $f\left(y_{ij},y_{ik}|\theta \right)$ by the univariate density functions of $y_{ij}$ and the difference $z_{ijk}=y_{ij}-y_{ik}.$ For the mean model (2.2), we have $y_{ij}~N\left(\mu ,{\sigma }_v^2+{\sigma }_e^2\right)$ and $z_{ijk}~N\left(0,2{\sigma }_e^2\right)$ . By reparametrizing $\theta ={\left(\mu ,{\sigma }_v^2,{\sigma }_e^2\right)}^T$ as $\text{ϕ}={\left(\mu ,{\sigma }^2,{\sigma }_e^2\right)}^T$ where ${\sigma }^2={\sigma }_v^2+{\sigma }_e^2,$ we see that the parameters of the two univariate density functions are distinct and the log composite likelihoods corresponding to $y_{ij}$ and $z_{ijk}$ are given by

$$l_{wCy}\left(\mu ,{\sigma }^2\right)={\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j\in s(i)}{w}_{j|i}\log f\left(y_{ij}|\mu ,{\sigma }^2\right)}}$$

and

$$l_{wCz}\left({\sigma }_e^2\right)={\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}\log f\left(z_{ijk}|{\sigma }_e^2\right).}}$$

We then solve the resulting weighted composite score equations

$${\widehat{U}}_{wCy1}\left(\mu ,{\sigma }^2\right)=\partial l_{wCy}\left(\mu ,{\sigma }^2\right)/\partial \mu ={\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j\in s_i}{w}_{j|i}\left(y_{ij}-\mu \right)/{\sigma }^2=0,}}$$

$${\widehat{U}}_{wCy2}\left(\mu ,{\sigma }^2\right)=\partial l_{wCy}\left(\mu ,{\sigma }^2\right)/\partial {\sigma }^2=\frac{1}{2}{\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j\in s(i)}{w}_{j|i}\left[-\frac{1}{{\sigma }_{}^2}+\frac{{\left(y_{ij}-\mu \right)}^2}{{\sigma }^4}\right]=0}}$$

$${\widehat{U}}_{wCz}\left({\sigma }_e^2\right)=\partial l_{wCz}\left({\sigma }_e^2\right)/\partial {\sigma }_e^2=\frac{1}{2}{\displaystyle \sum _i{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}\left(-\frac{1}{{\sigma }_e^2}+\frac{z_{ijk}^2}{2{\sigma }_e^4}\right)=0}}$$

to get the weighted composite likelihood (WCL) estimators ${\widehat{\mu }}_{wC},{\widehat{\sigma }}_{vwC}^2$ and ${\widehat{\sigma }}_{ewC}^2$ . The WCL estimators are identical to (3.9) $–$ (3.11) obtained by the weighted estimating equations approach of Section 3.

We now turn to the nested error linear regression model (3.13). We first note that $y_{ij}~N\left(x_{ij}^T\beta ,{\sigma }^2\right)$ where ${\sigma }^2={\sigma }_v^2+{\sigma }_e^2,$ and $z_{ijk}=y_{ij}-y_{ik}~N\left\{{\left(x_{ij}-x_{ik}\right)}^T\beta ,2{\sigma }_e^2\right\}.$ It follows that the weighted composite score equations are given by

$$\begin{array}{c}{\widehat{U}}_{wCy1}(\beta ,{\sigma }^2)=\partial l_{wCy}\left(\beta ,{\sigma }^2\right)/\partial \beta \\ ={\displaystyle \sum _{i\in s}{w}_i}{\displaystyle \sum _{j\in s(i)}{w}_{j|i}{x}_{ij}\left(y_{ij}-x_{ij}^T\beta \right)=0}\end{array}$$

$$\begin{array}{c}{\widehat{U}}_{wCy2}(\beta ,{\sigma }^2)=\partial l_{wCy}\left(\beta ,{\sigma }^2\right)/\partial {\sigma }^2\\ =-\frac{1}{2}{\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j\in s(i)}{w}_{j|i}\left[\frac{1}{{\sigma }^2}-\frac{{\left(y_{ij}-x_{ij}^T\beta \right)}^2}{{\sigma }^4}\right]=0}}\end{array}$$

and

$$\begin{array}{c}{\widehat{U}}_{wCz}({\sigma }_e^2)=\partial l_{wCz}\left({\sigma }_e^2\right)/\partial {\sigma }_e^2\\ =-\frac{1}{2}{\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}\left\{\frac{1}{{\sigma }_e^2}-\frac{{\left[z_{ijk}-{\left(x_{ij}-x_{ik}\right)}^T\beta \right]}^2}{2{\sigma }_e^4}\right\}=0.}}\end{array}$$

The resulting WCL estimators of $\beta$ , ${\sigma }_v^2$ and ${\sigma }_e^2$ are given by

$${\widehat{\beta }}_{wC}={\left({\displaystyle \sum _{i\in s}{\displaystyle \sum _{j\in s(i)}{w}_{ij}{x}_{ij}{x}_{ij}^T}}\right)}^{-1}\left({\displaystyle \sum _{i\in s}{\displaystyle \sum _{j\in s(i)}{w}_{ij}{x}_{ij}{y}_{ij}}}\right),$$

$${\widehat{\sigma }}_{wC}^2={\displaystyle \sum _{i\in s}{\displaystyle \sum _{j\in s(i)}{w}_{ij}{\left(y_{ij}-x_{ij}^T{\widehat{\beta }}_{wC}\right)}^2}}/{\displaystyle \sum _{i\in s}{\displaystyle \sum _{j\in s(i)}{w}_{ij}}},$$

and

$${\widehat{\sigma }}_{ewC}^2={\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}}}{\left[z_{ijk}-{\left(x_{ij}-x_{ik}\right)}^T{\widehat{\beta }}_{wC}\right]}^2/\left(2{\displaystyle \sum _{i\in s}{w}_i}{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}}\right).$$

The estimator of ${\sigma }_v^2$ is given by ${\widehat{\sigma }}_{vwC}^2={\widehat{\sigma }}_{wC}^2-{\widehat{\sigma }}_{ewC}^2.$ Again, the WCL estimators ${\widehat{\beta }}_{WC}$ , ${\widehat{\sigma }}_{vWC}^2$ and ${\widehat{\sigma }}_{eWC}^2$ are identical to (3.17)-(3.19) obtained from the weighted estimating equations approach of Section 3.

The above composite likelihood approach, based on $y_{ij}$ and $z_{ijk}=y_{ij}-y_{ik}$ , is not applicable to the linear two-level model given by (2.4) because the parameter vector, $\theta$ , is not identifiable under the composite likelihood obtained from the $y_{ij}$ and $z_{ijk}$ . We need the pairwise method to handle model (2.4).

Marginally, ${\left(y_{ij},y_{ik}\right)}^T$ is bivariate normal with means $x_{ij}^T\beta$ and $x_{ik}^T\beta$ and $2\times 2$ covariance matrix

$${\Sigma }_{i(jk)}=\left[\begin{array}{cc}{\sigma }_e^2+x_{ij}^T{\Sigma }_v{x}_{ij}& x_{ij}^T{\Sigma }_v{x}_{ik}\\ x_{ik}^T{\Sigma }_v{x}_{ij}& {\sigma }_e^2+x_{ik}^T{\Sigma }_v{x}_{ik}\end{array}\right].$$

It now follows from (4.3) that the weighted composite score equations are given by

$$\beta :{\widehat{U}}_{wC\beta }={\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}{X}_{i(jk)}^T{\Sigma }_{i(jk)}^{-1}\left(y_{i(jk)}-X_{i(jk)}^T\beta \right)}}=0\left(4.4\right)$$

and

$$\begin{array}{l}\tau :{\widehat{U}}_{wCl}=\frac{1}{2}{\displaystyle \sum _{i\in s}{w}_i{\displaystyle \sum _{j<k\in s(i)}{w}_{jk|i}[{\left(y_{i(jk)}-X_{i(jk)}^T\beta \right)}^T{\Sigma }_{i(jk)}^{-1}\frac{\partial {\Sigma }_{i(jk)}}{\partial {\tau }_l}{\Sigma }_{i(jk)}^{-1}\left(y_{i(jk)}-X_{i(jk)}^T\beta \right)}}\left(4.5\right)\\ -\text{tr}\left({\Sigma }_{i(jk)}^{-1}\frac{\partial {\Sigma }_{i(jk)}}{\partial {\tau }_l}\right)]=0,l=1,\mathrm{...},p(p+1)/2+1=P\end{array}$$

where $X_{i(jk)}$ is the $2\times p$ matrix with rows $x_{ij}^T$ and $x_{ik}^T$ , $y_{i(jk)}={\left(y_{ij},y_{ik}\right)}^T$ and $\tau$ is the P-vector with elements ${\tau }_1={\sigma }_e^2$ and the $p\left(p+1\right)/2$ distinct elements of ${\Sigma }_v$ denoted by ${\tau }_2,\mathrm{...},{\tau }_P$ . We can solve the weighted composite score equations (4.4) and (4.5) iteratively using the Newton-Raphson method or some other iterative method to obtain the WCL estimators ${\widehat{\beta }}_{wC}$ and ${\widehat{\tau }}_{wC}$ .

In the special case of the nested error linear regression model (3.13), the census score equations, based on the full census log-likelihood $l(\theta )$ given by (2.5), can be written in a closed form. The corresponding sample weighted score equations depend only on the level 1 weights $w_{j|i}$ and $w_{jk|i}$ and the level 2 weights $w_i,$ similar to the weighted composite score equations (see the Appendix). The resulting estimators are design-model consistent for $\theta$ , unlike the estimators based on the weighted pseudo log-likelihood $l_w(\theta )$ given by (2.7) and (2.8). However, for more complex models, such as two level models with random slopes, the sample weighted score equations will depend on third order and fourth order level 1 inclusion probabilities, unlike the weighted composite score equations (4.3) that depend only on the first order and second order level 1 inclusion probabilities, even for complex multi-level models. We have therefore not included the weighted score equations approach, based on the full census log-likelihood, in the simulation study.

Previous | Next