Suppose you were asked shortly after his inauguration in 2021, "Do you think Joe Biden will be effective as President of the United States?" Then, after four years, suppose you are asked, "Do you think Joe Biden has been effective as President of the United States?"

Whatever your responses are over the two time points, they will not be independent because they come from the same person. This is not the same situation as being chosen by chance in two random samples at those time points. If the same people are necessarily asked again four years later, then we really have only one sample with pairs of (dependent) responses.

In dependent samples, each observation in one sample can be paired with an observation in the other sample. Examples may include responses for essentially the same question at two points in time, such as in the example above, or even responses for different individuals if they are related in some way. Consider the examples below, the first of which has responses paired by their sibling relationship.

Let us now look at the example involving the siblings in more detail. Both siblings solve the same puzzle, and the response for each is whether it took less than or longer than one minute. It is sensible to assume that these two responses should be related because siblings likely inherit similar problem-solving skills from their common parents, and indeed if we test for independence between siblings, we have \(X^2=4.3612\) with one degree of freedom and p-value 0.03677. If we view solving the puzzle in less than one minute as "success", then this is equivalent to testing for equal row proportions.

However, the question of primary interest in this study is whether the older siblings tend to have a  higher proportion of success, compared with that of the younger siblings, which is a comparison of the first-row proportion against the first column proportion. Such a test does not require the use of siblings, and two samples of six-year-olds and eight-year-olds could have been independently chosen for this purpose. But using siblings allows for matched pairs of responses and controls for confounding factors that may be introduced with children from different parents.

The estimate for the difference in success proportions between ages is

\((15+7)/37-(15+5)/37=0.5946-0.5405=0.0541\)

Another way of asking the question of no age effect is to ask whether the margins of the table are the same (rows versus columns), which can be done with the test of marginal homogeneity or McNemar's test, which we look at next.

The notation needed for this test is the same as what we've seen earlier for the test of independence, but instead of focusing on comparing row proportions, we compare row versus column.

For older siblings, the probability of solving the puzzle in less than one minute (success) is

\(\pi_{1+} = \pi_{11} + \pi_{12}\)

And for younger siblings, the probability of solving the puzzle in less than one minute is

\(\pi_{+1} = \pi_{11} + \pi_{21}\)

The null hypothesis of no difference (marginal homogeneity) in a \(2 \times 2\) table is

\(H_0 \colon \pi_{1+} = \pi_{+1} \)

and is equivalent to the hypothesis that the off-diagonal probabilities are equal:

\(H_0 \colon \pi_{12} = \pi_{21} \)

The second of these above is also known as the hypothesis of symmetry and generalizes to equal off-diagonal cell counts for larger tables as well. Note that the diagonal elements are not important here. They correspond to the proportions of puzzles equally between the two siblings (either both less than or both greater than one minute). All the information required to determine whether there's a difference due to age is contained in the off-diagonal elements.

For general square \(I \times I \) tables, the hypothesis of marginal homogeneity is different from the hypothesis of symmetry, and the latter is a stronger hypothesis; symmetry introduces more structure in the square table. In a \(2 \times 2\) table, however, these two are the same test.

This is the usual test of marginal homogeneity (and symmetry) for a \(2 \times 2\) table.

\(H_0 : \pi_{1+} = \pi_{+1}\) or equivalently \(\pi_{12} = \pi_{21}\)

Suppose that we treat the total number of observations in the off-diagonal as fixed:

\(n^\ast =n_{12}+n_{21}\)

Under the null hypothesis above, each of \(n_{12}\) and \(n_{21}\) is assumed to follow a \(Bin (n^\ast , 0.5)\) distribution. The rationale is that, under the null hypothesis, we have \(n_{12}+n_{21}\) total "trials" that can either result in cell \((1,2)\) or \((2,1)\) with probability 0.5. And provided that \(n^*\) is sufficiently large, we can use the usual normal approximation to the binomial:

\(z=\dfrac{n_{12}-0.5n^\ast}{\sqrt{0.5(1-0.5)n^\ast}}=\dfrac{n_{12}-n_{21}}{\sqrt{n_{12}+n_{21}}}\)

where \(0.5n^*\) and \(0.5(1-0.5)n^\ast\) are the expected count and variance for \(n_{12}\) under the \(H_0\). Under \(H_0\), \(z\) is approximately standard normal. This approximation works well provided that \(n^* \ge 10\). The p-value would depend on the alternative hypothesis. If \(H_a\) is that older siblings have a greater success probability, then the p-value would be

\(P(Z\ge z)\)

where \(Z\) is standard normal, and \(z\) is the observed value of the test statistic. A lower-tailed alternative would correspondingly use the lesser-than probability, and a two-sided alternative would double the one-sided probability. Alternatively, for a two-sided alternative we may compare

\(z^2=\dfrac{(n_{12}-n_{21})^2}{n_{12}+n_{21}}\)

to a chi-square distribution with one degree of freedom. This test is valid under general multinomial sampling when \(n^\ast\) is not fixed, but the grand total \(n\) is. When the sample size is small, we can compute exact probabilities (p-values) using the binomial probability distribution.

Applying this to our example data gives

\(z=\dfrac{7-5}{\sqrt{7+5}}=0.577\)

The p-value is \(P(Z\ge 0.577)=0.2820\), which is not evidence of a difference in success probabilities (solving the puzzle in less than one minute) between the two age groups.

A sensible effect-size measure associated with McNemar’s test is the difference between the marginal proportions,

\(d=\pi_{1+}-\pi_{+1}=\pi_{12}-\pi_{21}\)

In large samples, the estimate of \(d\),

\(\hat{d}=\dfrac{n_{12}}{n}-\dfrac{n_{21}}{n}\)

is unbiased and approximately normal with variance

\begin{align}
V(\hat{d}) &= n^{-2} V(n_{12}-n_{21})\\
&= n^{-2}[V(n_{12})+ V(n_{21})-2Cov(n_{12},n_{21})]\\
&= n^{-1} [\pi_{12}(1-\pi_{12})+\pi_{21}(1-\pi_{21})+2\pi_{12} \pi_{21}]\\
\end{align}

An estimate of the variance is

\(\hat{V}(\hat{d})=n^{-1}\left[\dfrac{n_{12}}{n}(1-\dfrac{n_{12}}{n})+\dfrac{n_{21}}{n}(1-\dfrac{n_{21}}{n})+2\dfrac{n_{12}n_{21}}{n^2}\right]\)

and an approximate 95% confidence interval is

\(\hat{d}\pm 1.96\sqrt{\hat{V}(\hat{d})}\)

In our example, we get an estimated effect of \(\hat{d} = 0.0541\) and its standard error of \(\sqrt{\hat{V}(\hat{d})}=0.0932\), giving 95% confidence interval

\(0.0541\pm 1.96(0.0932)=(-0.1286,  0.2368)\)

Thus, although the older siblings had a higher proportion of success, it was not statistically significant. We cannot conclude that the two-year age difference is associated with faster puzzle-solving times.

Next, we do this analysis in SAS and R.

For the sibling puzzle-solving example, note that the data consisted of 37 responses for six-year-olds (younger siblings) and 37 responses for eight-year-olds (older siblings). How would the results differ if, instead of siblings, those same values had arisen from two independent samples of children of those ages?

To see why the approach with the dependent data, matched by siblings, is more powerful, consider the table below. The sample size here is defined as \(n = 37+37=74\) total responses, compared with \(n=37\) total pairs in the previous approach.

The estimated difference in proportions from this table is identical to the previous one:

\(\hat{d}=\hat{p}_1-\hat{p}_2=\dfrac{22}{37}-\dfrac{20}{37}=0.0541\)

But with independent data, the test for an age effect, which is equivalent to the usual \(\chi^2\) test of independence, gives \(X^2 = 0.2202\), compared with McNemar's test gave \(z^2 = 0.333\). The standard error for \(\hat{d}\) in this new table is

\(\sqrt{\dfrac{\hat{p}_1(1-\hat{p}_1)}{37}+\dfrac{\hat{p}_2(1-\hat{p}_2)}{37}}=0.1150\)

compared with \(0.0932\), when using siblings and taking into account the covariance between them. Just as with matched pairs for quantitative variables, the covariance between pair members leads to smaller standard errors and greater overall power, compared with the independent samples approach.

data matched;
input time $ approval $ count ;
datalines;
 t1 approve 944
 t1 disapprove 656
 t2 approve  880
 t2 disapprove 720
;
proc freq data=matched order=data;
 weight count;
    tables time*approval /chisq riskdiff;
run;

notsiblings = matrix(c(22,20,15,17),nr=2,
 dimnames=list(c("Older","Younger"),c("<1 min",">1 min")))
notsiblings

chisq.test(notsiblings, correct=F)
prop.test(notsiblings, correct=F)

> chisq.test(notsiblings, correct=F)

        Pearson's Chi-squared test

data:  notsiblings
X-squared = 0.22024, df = 1, p-value = 0.6389

> prop.test(notsiblings, correct=F)

        2-sample test for equality of proportions without continuity
        correction

data:  notsiblings
X-squared = 0.22024, df = 1, p-value = 0.6389
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.1713610  0.2794691

McNemar’s test applies whenever the hypothesis of marginal homogeneity is of interest. Dependency and marginal homogeneity may arise in varieties of problems involving dependency.

Another hypothesis of interest is to evaluate to what degree two different examiners agree on different systems of evaluation. This has important applications in medicine where two physicians may be called upon to evaluate the same group of patients for further treatment.

The Cohen's Kappa statistic (or simply kappa) is intended to measure agreement between two variables.

A rather restrictive form of agreement between two variables is symmetry, where the expected counts satisfy \(\mu_{ij}=\mu_{ji}\) for all \(i\) and \(j\). As the name suggests, this means the table of expected counts is a symmetric matrix.

Using the notation from the Siskel and Ebert example with general log-linear model

\(\log\mu_{ij}=\lambda+\lambda^S_i+\lambda^E_j+\lambda^{SE}_{ij}\)

the symmetry assumption would mean that

• \(\lambda_{ij}^{SE}=\lambda^{SE}_{ji}\)
• \(\lambda^S_i=\lambda^E_i\)

for all \(i\) and \(j\). Note that this differs from the independence assumption that all \(\lambda_{ij}^{SE}=0\). Both are restrictive models, but neither is a special case of the other.

For an \(I\times I\) table, then there are \(I\) total counts along the main diagonal and \(I(I-1)\) total counts off the main diagonal, giving \(I+(I(I-1)/2\) unique parameters to be estimated to fit the symmetry model.

The challenge in actually fitting this model, however, is that the symmetry restrictions are not simply removing certain terms (setting them equal to 0). Rather, the restriction is that certain counts---the off-diagonal ones---are required to be equal. So, we use a regression approach with some carefully defined indicator variables. Specifically, for each observation in the sample, we take \(I_{ij}=1\) if the observation falls in row \(i\) and column \(j\) OR if the observation falls in row \(j\) and column \(i\), for all \(i\le j\le I\); otherwise \(I_{ij}=0\).

The idea is that the symmetry model views the \(n_{ij}\) count and the \(n_{ji}\) count as coming from a common population, and so only a single parameter is attributed to its mean. For the Siskel and Ebert example, this reduces to

\(\log\mu_{ij}=\beta_1I_{11}+\beta_2I_{22}+\beta_3I_{33}+ \beta_4I_{12}+\beta_5I_{13}+\beta_6I_{23} \)

Thus, the \(i^{th}\) main diagonal mean is \(\mu_{ii}=e^{\beta_i}\), and the \((i,j)^{th}\) off-diagonal mean is \(\mu_{ij}=\mu_{ji}=e^{\beta}\), where \(\beta\) is the coefficient for \(I_{ij}\) (for \(i<j\)).

Let's see how we can do this in SAS and R.

data critic;
set critic;

/* indicator variables for each diagonal cell*/
I11=0; if siskel='con' AND ebert='con' THEN I11=1;
I22=0; if siskel='mixed' AND ebert='mixed' THEN I22=1;
I33=0; if siskel='pro' AND ebert='pro' THEN I33=1;

/* indicator variables for each off-diagonal cell*/
I12=0; if siskel='con' AND ebert='mixed' THEN I12=1;
 if siskel='mixed' AND ebert='con' THEN I12=1;
I13=0; if siskel='con' AND ebert='pro' THEN I13=1;
 if siskel='pro' AND ebert='con' THEN I13=1;
I23=0; if siskel='mixed' AND ebert='pro' THEN I23=1;
 if siskel='pro' AND ebert='mixed' THEN I23=1;
; run;

/* symmetry */
proc genmod data=critic order=data;
model count=I11 I22 I33 I12 I13 I23 / 
 link=log dist=poisson noint predicted;
title "Symmetry Model";
run;

S = factor(rep(c("con","mixed","pro"),each=3)) # siskel
E = factor(rep(c("con","mixed","pro"),times=3)) # ebert
count = c(24,8,13,8,13,11,10,9,64)

# indicator variables for diagonals
I11 = (S=="con")*(E=="con")
I22 = (S=="mixed")*(E=="mixed")
I33 = (S=="pro")*(E=="pro")

# indicator variables for off-diagonals
I12 = (S=="con")*(E=="mixed") + (E=="con")*(S=="mixed")
I13 = (S=="con")*(E=="pro") + (E=="con")*(S=="pro")
I23 = (S=="mixed")*(E=="pro") + (E=="mixed")*(S=="pro")

# symmetry model
fit.s = glm(count~I11+I12+I13+I22+I23+I33-1,family=poisson(link='log'))

> summary(fit.s)

Call:
glm(formula = count ~ I11 + I12 + I13 + I22 + I23 + I33 - 1, 
    family = poisson(link = "log"))

Deviance Residuals: 
      1        2        3        4        5        6        7        8  
 0.0000   0.0000   0.4332   0.0000   0.0000   0.3112  -0.4525  -0.3217  
      9  
 0.0000  

Coefficients:
    Estimate Std. Error z value Pr(>|z|)    
I11   3.1781     0.2041  15.569   <2e-16 ***
I12   2.0794     0.2500   8.318   <2e-16 ***
I13   2.4423     0.2085  11.713   <2e-16 ***
I22   2.5649     0.2774   9.248   <2e-16 ***
I23   2.3026     0.2236  10.297   <2e-16 ***
I33   4.1589     0.1250  33.271   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 721.15855  on 9  degrees of freedom
Residual deviance:   0.59276  on 3  degrees of freedom
AIC: 52.784

Number of Fisher Scoring iterations: 4

> mu.s = fitted.values(fit.s)
> matrix(mu.s,nrow=3)
     [,1] [,2] [,3]
[1,] 24.0    8 11.5
[2,]  8.0   13 10.0
[3,] 11.5   10 64.0

Although the symmetry model fits well in this example, it may be too restrictive for other types of data. So, something less restrictive (requiring less agreement) may be more suitable. This leads us to consider "quasi-symmetry".

Less restrictive than the symmetry model but still assuming some agreement compared with the fully saturated (unspecified) model, the quasi-symmetry model assumes that

\(\lambda_{ij}^{SE}=\lambda^{SE}_{ji}\)

Note that, compared with the symmetry model, the quasi-symmetry model does not require marginal homogeneity, which is that \(\lambda^S_i=\lambda^E_i\). In terms of the notation we have so far, the model now becomes

\(\log(\mu_{ij}) = \lambda+\lambda_i^S+\lambda_j^E+ \beta_4I_{12}+\beta_5I_{13}+\beta_6I_{23} \)

Thus, we use the separate S and E factors to allow for more general marginal distributions but apply a restriction to their joint distribution.

Generally for an \(I\times I\) table, the number of parameters breaks down as follows:

\(1+(I-1)+(I-1)+\dfrac{I(I-1)}{2}\)

The first three quantities allow one parameter for \(\lambda\) and \(I-1\) parameters for each marginal distribution separately for S and E. The last quantity \(I(I-1)/2\) corresponds to the number of off-diagonal counts and is identical to the symmetry model. But the symmetry model requires only \(I\) parameters for a common marginal distribution, whereas the quasi-symmetry model here allows for additional parameters and is hence less restrictive.

While it may seem at a glance that using the same indicators for the off-diagonal table counts would require symmetry, note that the additional parameters allowed in the marginal distributions for S and E influence the off-diagonal counts as well. For example, the expected count for the \((1,2)\) cell would be

\(\mu_{12}=\exp[\lambda+\lambda_1^S+\lambda_2^E+\beta_4]\)

while the expected count for the \((2,1)\) cell would be

\(\mu_{21}=\exp[\lambda+\lambda_2^S+\lambda_1^E+\beta_4]\)

For the \(3\times3\) table, we're currently working with, the quasi-symmetry model has 8 parameters, and so the deviance test versus the saturated model will have only a single degree of freedom. But for larger tables, it can be a more moderate appreciable reduction. As it is, we see below when fitting this model with software, the test statistic is \(G^2=0.0061\) and is very similar to the saturated model.

/* quasi-symmetry */
proc genmod data=critic order=data;
class siskel ebert;
model count=siskel ebert I12 I13 I23 /
 link=log dist=poisson predicted;
title "Quasi-Symmetry Model";
run;

# quasi-symmetry model
fit.q = glm(count~S+E+I12+I13+I23,family=poisson(link='log'))

> summary(fit.q)

Call:
glm(formula = count ~ S + E + I12 + I13 + I23, family = poisson(link = "log"))

Deviance Residuals: 
       1         2         3         4         5         6         7         8  
 0.00000  -0.03454   0.02727   0.03482   0.00000  -0.02949  -0.03092   0.03281  
       9  
 0.00000  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.1781     0.2041  15.569  < 2e-16 ***
Smixed       -0.3188     0.2594  -1.229  0.21913    
Spro          0.3679     0.2146   1.714  0.08652 .  
Emixed       -0.2943     0.2594  -1.134  0.25666    
Epro          0.6129     0.2146   2.856  0.00430 ** 
I12          -0.7921     0.3036  -2.609  0.00907 ** 
I13          -1.2336     0.2414  -5.110 3.22e-07 ***
I23          -1.0654     0.2712  -3.928 8.55e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 1.0221e+02  on 8  degrees of freedom
Residual deviance: 6.0515e-03  on 1  degrees of freedom
AIC: 56.198

> mu.q = fitted.values(fit.q)
> matrix(mu.q,nrow=3)
          [,1]      [,2]      [,3]
[1,] 24.000000  7.901913 10.098087
[2,]  8.098087 13.000000  8.901913
[3,] 12.901913 11.098087 64.000000

We have seen this when we discussed the two-way table for the siblings' puzzle times. Marginal homogeneity means the row and column distributions (of a square table) are the same. That is, 

\(\mu_{i+}=\mu_{+i}\)

for all \(i\). There is no log-linear model that corresponds directly to marginal homogeneity, but we can use the log-linear models of symmetry and quasi-symmetry to indirectly test for marginal homogeneity.

Recall that symmetry has two restrictions:

• \(\lambda_{ij}^{SE}=\lambda^{SE}_{ji}\)
• \(\lambda^S_i=\lambda^E_i\)

And this suggests that symmetry consists of restrictions from both quasi-symmetry and marginal homogeneity. This correspondence is actually one-to-one, meaning that both quasi-symmetry and marginal homogeneity together are necessary and sufficient for symmetry. And thus if a symmetry model lacks fit when compared with a quasi-symmetry model, it must be due to the marginal homogeneity assumption. This gives us an avenue for testing marginal homogeneity:

\(G^2\) (marginal homogeneity versus saturated) \(= G^2\) (symmetry versus quasi-symmetry) \)

For our example, we calculated \(G^2 = 0.5928\) for the test of symmetry versus saturated (with 3 df) and \(G^2=0.0061\) for the test of quasi-symmetry saturated (with 1 df), giving a likelihood ratio test of \(G^2=0.5928-0.0061=0.5867\) for the test of symmetry versus quasi-symmetry (with 2 df). Indirectly, this serves as our test for marginal homogeneity, and in this case, it would be a reasonable fit to the data (p-value 0.746).

A diagram of the models discussed in the lesson are depicted below.

The most general (complex) model is the saturated model. Each of quasi-symmetry and marginal homogeneity is a special case of the saturated model, but neither is a special case of the other. And symmetry is both quasi-symmetry and marginal homogeneity.

In terms of deviance statistics versus saturated, the table below summarizes the test results. For this example, the symmetry model would be preferred because it's the most restrictive model that fits the data reasonably well.

Thus, Siskel and Ebert tend to agree on their ratings of movies.