The idea behind GEEs is to produce reasonable estimates of model parameters, along with standard errors, without specifying a likelihood function in its entirety, which can be quite difficult with a multivariate categorical response. We have to account for the correlation among the multiple responses that arise from a single subject, but we can largely estimate these correlations from the data without having to accurately specify the form of the correlation structure.

In the General Social Survey (GSS) of 2018, respondents answered questions regarding their confidence in various institutions in the U.S., including education, medicine, and the scientific community. Among these three institutions, are there differences in the proportions of people having "a great deal" of confidence? Responses are clustered by subject because each subject provided answers to all three questions.

The response vector for subject \(i\) is \(Y_i=(Y_{i1},\ldots,Y_{in_i})\), where \(Y_{ij}\) is the categorical response for subject \(i\) (\(1,\dots,n\)) and measurement \(j\) (\(1,\ldots,n_i\)). With expected value \(E(Y_{i})=\mu_i'=(\mu_{i1},\ldots,\mu_{in_i})\) and covariate/predictor vector \(x_{ij}\), we have the general model expression:

\(g(\mu_{ij})=x_{ij}'\beta\)

In our example above, \(Y_{ij}\) is binomial with mean \(\mu_{ij}=\pi_{ij}\), and the logit link would be used for \(g\). If the institution indicators, say \(Med_{ij}=1\) for medicine and \(Sci_{ij}=1\) for scientific community (each would be 0 otherwise), are added, along with a covariate for age, we have \(x_{ij}'=(1,Med_{ij},Sci_{ij},Age_i)\) and

\(\log\left(\dfrac{\pi_{ij}}{1-\pi_{ij}}\right)=\beta_0+Med_{ij}\beta_1+Sci_{ij}\beta_2+Age_i\beta_3\)

Note that the age of the subject does not vary within subject and so does not require a \(j\) subscript. Although subjects are assumed to be independently collected, the multiple responses within a subject are allowed to have a correlation structure. Here are some potential structures we may consider. In each matrix below (in general, these would be \(n_i\times n_i\)), the value in the \(j^{th}\) row and \(k^{th}\) column represents the correlation between the \(j^{th}\) and \(k^{th}\) measurement for the same subject.

Maximum-likelihood estimation of these parameters from the would be identical to that for the familiar logistic model, except that the correlation structure would be incorporated into the likelihood for responses from the same subject. The GEE alternative avoids this.

Interpretation of Parameter Estimates:

Model Fit:

Advantages:

Limitations:

Starting with \(E(y_{i})=\mu_{i}\), the vector of means for subject \(i\) connected with the predictors via \(g(\mu_i)=x_i'\beta)\), we let \(\Delta_i\) be the diagonal matrix of variances

\(\Delta_i=\text{Diag}[\text{Var}(y_{ij})]= \begin{bmatrix}
Var_{i1} &\cdots & \cdots &\vdots \\
\vdots & Var_{i2} & \cdots & \vdots\\
\vdots& \cdots & \ddots & \vdots\\
\vdots& \cdots & \cdots & Var_{ij}
\end{bmatrix}\).

In terms of the correlation structure, this would be:

\( \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\)

Unless the observations within a subject are independent,

\(\Delta_i \neq \text{Cov}(y_i)\)

But if \(\Delta_i\) were correct, we could stack up the \(y_i\)'s and estimate \(\beta\) using techniques for generalized linear models.

How bad is it to pretend that \(\Delta_i\) is correct? Let \(\hat{\beta}\) be the estimate that assumes observations within a subject are independent (e.g., as found in ordinary linear regression, logistic regression, etc.)

• If \(\Delta_i\) is correct, then \(\hat{\beta}\) is asymptotically unbiased and efficient. Recall that unbiased \(E(\hat{\beta})=\beta\), efficient means it has the smallest variance of all other possible estimations.
• If Δ_i is not correct, then \(\hat{\beta}\) is still asymptotically unbiased but no longer efficient

  \(\hat{\sigma}^2 \left[X^T \hat{W} X \right]^{-1}\)

  may be very misleading (here, X is the matrix of stacked X_i's and \(\hat{W}\) is the diagonal matrix of final weights, if any)

  • The 'naive' standard error for \(\hat{\beta}\), obtained from the naive estimate of \(\text{Cov}(\hat{\beta})\)
  • consistent standard errors for \(\hat{\beta}\) are still possible using the sandwich estimator (sometimes called the 'robust' or 'empirical' estimator)

The sandwich estimator was first proposed by Huber (1967) and White(1980); Liang and Zeger (1986) applied it to longitudinal data

\(\left[X^T \hat{W} X \right]^{-1} \left[\sum\limits_i X_i^T (y_i-\hat{\mu}_i)(y_i-\hat{\mu}_i)^T X_i \right]\left[X^T \hat{W} X \right]^{-1}\)

• provides a good estimate of \(\text{Cov}(\hat{\beta})\) in large samples (several hundred subjects or more) regardless of the true form of \(\text{Cov}\left(y_i\right)\)
• in smaller samples it could be rather noisy, so 95% intervals obtained by ± 2 SE's could suffer from under-coverage

When within-subject correlations are not strong, Zeger (1988) suggests that the use of IEE with the sandwich estimator is highly efficient.

With PROC GENMOD, we can also try alternative assumptions about the within-subject correlation structure.

The type=exch or type=cs option specifies an "exchangeable" or "compound symmetry assumption," in which the observations within a subject are assumed to be equally correlated:

\(Corr(y_i)=\begin{bmatrix}
1 & \rho & \rho & \cdots & \rho\\
\rho & 1 & \rho & \cdots & \rho\\
\rho & \rho & 1 & \cdots & \rho\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
\rho & \rho & \rho & \cdots & 1
\end{bmatrix}\)

This is the same kind of structure that arises in a normal linear model with random intercepts. It acknowledges that correlations exist within a subject, but it assumes that the correlations are constant over time.

An autoregressive(AR-1) structure, specified by type=ar, allows the correlations to diminish over time as

\(Corr(y_{ij},y_{ij'})=\rho ^{|t_{ij}-t_{ij'}|}\),

where \(t_{ij}\) and \(t_{ij′}\) are the observation times for \(y_{ij}\) and \(y_{ij'}\) .

The exchangeable and the autoregressive structures both express the within-subject correlations in terms of a single parameter \(\rho\). If the subjects are measured at a relatively small common set of occasions, we may be able to estimate an arbitrary correlation matrix. With \(n_i = 4\) measurement times per subject, the unstructured matrix would have six correlations to estimate. An unstructured matrix is obtained by the option type=un.

 Which correlation structure is best for any given problem?

In principle, it makes sense to think that the one that is most nearly "correct" would be best. For example, if the autoregressive or exchangeable assumptions are wrong, one might think that the unstructured form would provide the most efficient answer. But that is not necessarily true, because moving to an unstructured matrix introduces many more unknown parameters which could destabilize the model fit. In certain cases, a wrong structure with a small number of parameters could perform better than a correct structure with many parameters. Unfortunately, the literature does not provide much guidance on how to choose a good working correlation structure for GEE.

Suppose that we observe responses \(Y_1,Y_1 , \dots , Y_N\) which we want to relate to covariates. The ordinary linear regression model is

\(Y_i \sim N(x_i^T\beta,\sigma^2)\)

where \(x_i\) is a vector of covariates, \(\beta\) is a vector of coefficients to be estimated, and \(\sigma^2\) is the error variance. Now let's generalize this model in two ways:

Introduce a link function for the mean \(E ( Y_i ) = \mu_i\),

\(g(\mu_i)=x_i^T\beta\).

We could also write this as

\(E(Y_i)=\mu_i(\beta)\quad\) (2)

where the covariates and the link become part of the function \(\mu_i\). For example, in a loglinear model we would have \(\mu_i=\exp(x_i^T\beta)\).

Allow for heteroscedasticity, so that

\(\text{Var}(Y_i)=V_i\quad\) (3) 

In many cases \(V_i\) will depend on the mean \(\mu_i\) and therefore on \(\beta\), and possibly on additional unknown parameters(e.g. a scale factor). For example, in a traditional overdispersed loglinear model, we would have \(V_i=\sigma^2 \mu_i=\sigma^2 \exp(x_i^T\beta)\).

A maximum-likelihood estimate for \(\beta\) under this model could be computed by a Fisher scoring procedure. ML estimates have two nice theoretical properties: they are approximately unbiased and highly efficient. Interestingly, the asymptotic theory underlying these properties does not really depend on the normality of \(Y_i\), but only on the first two moments. That is if the mean function (2) and the variance function (3) are correct, but the distribution of \(Y_i\) is not normal, the estimate of \(\beta\) obtained by maximizing the normal loglikelihood from

\(Y_i \sim N(\mu_i(\beta), V_i(\beta))\)

is still asymptotically unbiased and efficient.

If we maximize the normality-based loglikelihood without assuming that the response is normally distributed, the resulting estimate of \(\beta\) is called a quasi-likelihood estimate. The iterative procedure for computing the quasi-likelihood estimate, called quasi-scoring, proceeds as follows.

First, let's collect the responses and their means into vectors of length \(N\),

\(Y=
\begin{bmatrix}
Y_1\\
Y_2\\
\vdots\\
Y_N
\end{bmatrix},
\mu=
\begin{bmatrix}
\mu_1\\
\mu_2\\
\vdots\\
\mu_N
\end{bmatrix} \)

Also, let \(V\) be the \(N \times N\) matrix with \(V_1 , \dots , V_N\) on the diagonal and zeros elsewhere. Finally, let

\(D=\dfrac{\partial \mu}{\partial \beta}\)

be the \(N \times p\) matrix containing the partial derivatives of \(\mu\) with respect to \(\beta\). That is, the \(\left(i, j\right)^{th}\) element of \(D\) is \(\dfrac{∂\mu_i }{∂\beta_i}\) . In a simple linear model with an identity link, \(\mu_i=x_i^T\beta\), \(D\) is simply the \(N \times p\) matrix of covariates

\(X=
\begin{bmatrix}
x_1^T\\
x_2^T\\
\vdots\\
x_N^T
\end{bmatrix}\)

The iterative quasi-scoring procedure is

\(\beta_{\mbox{new}}-\beta_{\mbox{old}}= \left( D^T V^{-1} D\right)^{-1} D^T V^{-1}(Y-\mu)\quad\) (4)

where on the right-hand side of (4) the quantities \(D\), \(\tilde{V}\), and \(\mu\) are evaluated at \(\beta = \beta_{\mbox{old}}\).

It is easy to see that in the case of simple linear regression with homoscedastic response (\(\mu_i=x_i^T\beta\) and \(V_i=\sigma^2\) ), the quasi-scoring procedure converges to the OLS estimate

\(\hat{\beta}=\left(X^T X\right)^{-1}X^T Y\)

in a single step regardless of \(\beta_{\mbox{old}}\). In other cases (e.g. loglinear regression where \(\mu_i=\exp(x_i^T\beta)\), \(V_i=\sigma^2 \mu_i\) ) it produces the same estimate for \(\beta\) that we obtained earlier by fitting the generalized linear model.

Quasi-likelihood estimation is really the same thing as generalized linear modeling, except that we no longer have a full parametric model for \(Y_i\).

If we let \(U\) be the \(p \times 1\) vector

\(U=D^TV^{-1}(Y-\mu)\)

the final estimate from the quasi-scoring procedure satisfies the condition \(U = 0\). \(U\) can be regarded as the first derivative of the quasi-loglikelihood function. The first derivative of a loglikelihood is called a score, so \(U\) is called a quasi-score. \(U = 0\) is often called the set of estimating equations, and the final estimate for \(\beta\) is the solution to the estimating equations.

Estimating equations with a working variance function. We'll suppose that the mean regression function \(\mu_i \left(\beta\right)\) has been correctly specified but the variance function has not. That is, the data analyst incorrectly supposes that the variance function for \(Y_i\) is \(\tilde{V}_i\) rather than \(V_i\) , where \(\tilde{V}_i\) is another function of \(\beta\). The analyst then estimates \(\beta\) by solving the quasi-score equations

\(U=D^T \tilde{V}^{-1}(Y-\mu)=0\)

where \(\tilde{V}_i\) is the diagonal matrix of \(\tilde{V}_i\)s. Obviously it would be better to perform the scoring procedure using the true variance function \(V = \text{Diag}\left(V_1 , \dots , V_N\right)\) rather than \(\tilde{V}\). What are the properties of this procedure where \(\tilde{V}\) has been used instead of \(V\)?

1. The estimate \(\hat{\beta}\) is a consistent and asymptotically unbiased estimate of \(\beta\), even if \(\tilde{V} \neq V\). It's also asymptotically normal.
2. If \(\tilde{V} \neq V\) then \(\hat{\beta}\) is not efficient; that is, the asymptotic variance of \(\hat{\beta}\) is lowest when \(\tilde{V} = V\).
3. If \(\tilde{V}= V\) then the final value of the matrix \(\left(D^T V^{-1} D\right)^{-1}\) from the scoring procedure (4) (i.e. the value of this matrix with \(\hat{\beta}\) substituted for \(\beta\)) is a consistent estimate of \(\text{Var}(\hat{\beta})\).
4. If \(\tilde{V} \neq V\) then the final value of the matrix \(\left(D^{T} V^{-1} D\right)^{-1}\) from (4) is not a consistent estimate of \(\text{Var}(\hat{\beta})\).

Because of these properties, \(\hat{\beta}\) may still be a reasonable estimate of \(\beta\) if \(\tilde{V} \neq V\), but the final value of \((D^T \tilde{V}^{-1}D)^{-1}\) --- often called the "model-based" or "naive" estimator --- will not give accurate standard errors for the elements of \(\hat{\beta}\). However, this problem can be corrected by using the "robust" or "sandwich estimator," defined as

\(\left(D^T \tilde{V}^{-1}D\right)^{-1} \left(D^T \tilde{V}^{-1} E \tilde{V}^{-1} D\right) \left(D^T \tilde{V}^{-1}D\right)^{-1}\quad\), (5)

where,

\(E=\text{Diag}\left((Y_1-\mu_1)^2,\ldots, (Y_N-\mu_N)^2\right)\quad\) , (6)

and all quantities in (5) and (6) are evaluated with \(\hat{\beta}\) substituted for \(\beta\). The sandwich estimator is a consistent estimate of \(\text{Var}(\hat{\beta})\) even when \(\tilde{V} \neq V\).

The sandwich estimator was first proposed by Huber (1967) and White (1980), but was popularized in the late 1980s when Liang and Zeger extended it to multivariate or longitudinal responses.

As we consider higher-dimensional contingency tables, we are more likely to encounter sparseness. A sparse table is one where there are many cells with small counts and/or zeros. How many is 'many' and how small is 'small' are relative to

• the sample size \(n\)
• the table size \(N\)

We can have a small sample size, or a large sample size with a large number of cells. The large number of cells can be either due to the large number of classification variables, or small number of variables but with lots of levels.

Usually, when we talk about sparseness we are talking about zeros. This is a problem that is frequently encountered. There are different types of zeros, and we need to differentiate between the two.

Sampling Zeros occur when there is no observation in a cell, but there is a positive probability of observing an observation. That is, \(n_{ij} = 0\), but \(\pi_{ij} > 0\), for some \(i\) and \(j\). And increasing the sample size will eventually lead to an observation there.

Structural Zeros are are theoretically impossible to observe. That is, both \(n_{ij} = 0\), and \(\pi_{ij} = 0\). Tables with structural zeros are structurally incomplete. They are also known as incomplete tables. This is different from a partial classification where an incomplete table results from not being able to completely cross-classify all individuals or units.

Sample of adults on whether they got an infection once or twice:

The explanation here is that one cannot get an infection the second time until one gets it the first time. Thus, \(n_{21}\) is a structural zero.

If there is a single sampling zero in a table, the MLE estimates are always positive. If there are non-zero margins, but a bad pattern of zeros, some MLE estimates of the cells counts can be negative.

If there are zero margins, we cannot estimate the corresponding term in the model, and the partial odds ratios involving this cell will equal 0 or 1. For example, if \(n_{+11} = 0\), there is no MLE estimate for \(\lambda_{11}^{YZ} = 0\), and any model including this term will not fit properly. We can force \(\lambda_{11}^{YZ} = 0\) but then must adjust for the degrees of freedom in the model, and no software will do this automatically.

A general formula for adjusting the degrees of freedom:

\(df = \left(T_e - Z_e\right) - \left(T_p - Z_p\right)\)

where \(T_e \) is the number of cells fitted, \(Z_e\) is the number of zero cells, \(T_p\) is the number of parameters in the model, and \(Z_p \) is the number of parameters that cannot be estimated because of zero margins.

Let's see how this works in the example below. The data, also available in Zeros.sas and Zeros.R, are from Fienberg (1980) and contain a marginal total if we sum over X:

Z=1:

Z=2

The three marginal tables are

XY:

XZ:

YZ:

Here, \(n_{+11} = 0\), so there is no MLE estimate for \(\lambda_{11}^{YZ} = 0\), and any model including this term will not fit properly. Let's take a look at these examples in R and SAS.

There are two major effects on modeling contingency tables with sampling zeros:

• Maximum likelihood estimates (MLEs) may not exist for loglinear/logit models
• If MLE estimates exist, they may be biased

Recall that terms in a log-linear model correspond to the margins of the tables. Existence of MLE estimates depends on the positions of zeros in the table and what effects are included in a model.

If all \(n_{ij} > 0\) than MLE estimates of parameters are finite. If a table has a 0 marginal frequency, and there is a term in the model corresponding to that margin, the MLE estimates are infinite.

As mentioned above, the ML estimates of odds ratio can be severely biased. Also, the sampling distributions of fit statistics may be poorly approximated by the chi-squared distribution, which is why the standard errors are so large.

One suggested ad hoc solution is to add .5 to each cell in the table. Doing so shrinks the estimated odds ratios that are 1 to finite values and increases estimates that would otherwise be 0.

Other suggested solutions:

• When a model does not converge, try adding a tiny number to all zero cells in the table, such as 0.000000001. Some argue you should add this value to all of the cells. Essentially, adding very small values is what the software programs like SAS and R are doing when they fit these types of models.
• Extend this by doing a sensitivity analysis by adding different numbers of varying sizes to the cells. Examine fit statistics and parameter estimates to see if they change very much.

Unfortunately there is no single rule that covers all situations. We have to be careful about application of these situation-specific rules.

So while the \(G^2\) statistic for model fit may not be well approximated by the chi-squared distribution, the difference between \(G^2\)s for two nested models will be closer to chi-squared than the \(G^2\) for either model versus the saturated one. Chi-squared comparison tests depend more on the size of marginals than on cell sizes in the joint table. So if margins have cells exceeding 5, the chi-squared approximation of \(G^{2}\left(\mbox{reduced model}\right) - G^{2}\left(\mbox{full model}\right)\) should be reasonable.

• When \(df > 1\), it is OK to have the \(\hat{\mu}\) as small as 1 as long as fewer than 20% of the cells have \(\hat{\mu} < 5\)
• The permissible size of \(\hat{\mu}\) decreases as the size of the table \(N\) increases.
• The chi-squared distribution of \(X^{2}\) and \(G^{2}\) can be poor for sparse tables with both very small and very large \({\hat{\mu}}'s\) (relative to \(n/N\)).
• \(G^2\) is usually poorly approximated by the chi-squared distribution when \(\dfrac{n}{N} < 5\).The p-values for \(G^2\) may be too large or too small (it depends on n/N).
• For fixed \(n\) and \(N\), chi-squared approximations are better for smaller df than for larger df.

In this lesson, we saw an alternative way to estimate parameters and their standard errors for a GLM that doesn't rely on full specification of the likelihood function. This "generalized estimating equation" approach needs only the marginal distribution of the responses, along with a correlation structure to describe pairwise associations between observations within the same subject.

Additionally, we saw different ways in which 0s can appear in categorical data, depending on whether an observation could have been observed but wasn't or whether an observation could not have been observed at all due to its nature.