1 License

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (http://creativecommons.org/licenses/by-sa/4.0/).

The R markdown source for this document is https://www.stat.umn.edu/geyer/5421/notes/ch11.Rmd.

2 R

library(CatDataAnalysis)
library(ump)
library(alabama)
## Loading required package: numDeriv

3 Matched Pairs and Dichotomous Data

Up to now, when we had a contingency table, each count in each cell of the table was the classification of one individual, and the individuals were IID.

Now we are going to analyze data for which the individuals are still IID (so we are still in an IID situation), but for each individual we have two dichotomous responses (before and after, right and left, treatment and control, question 1 and question 3, whatever the two responses may be). Thus we get a two-by-two table.

Agresti’s example is the data in his Table 11.1

y <- matrix(c(175, 54, 16, 188), 2, 2)
dimnames(y) <- list(Election2004 = c("Democrat", "Republican"),
    Election2008 = c("Democrat", "Republican"))
y
##             Election2008
## Election2004 Democrat Republican
##   Democrat        175         16
##   Republican       54        188

These are data on 216.5 individuals on each of which we have two votes recorded. The vast majority of individuals voted the same way in both elections: 175 voted Democratic twice and 188 voted Republican twice. As we shall see, the other two cells of the table are the important ones. 16 individuals voted Democratic in 2004 and switched to Republican in 2008, and 54 switched the other way.

The responses within pairs are clearly dependent. As is usual with matched pairs we just assume this dependence rather than test for it. We could also say we are allowing for any dependence (which is expected), but our test is valid whether there is any dependence or not.

We can also use the procedures for matched pair dichotomous data even when the individuals are only “individuals” (in scare quotes). In so-called case-control studies individuals receiving treatment are matched with different individuals not receiving treatment who match the treated individuals on certain covariates (Agresti, Section 11.2.5).

4 McNemar’s Test

As the first displayed equation in Section 11.1 in Agresti says, when we compare population proportions before and after (or whatever the matched pairs may be) we have \[\begin{equation} \tag{$*$} p_{1 +} - p_{+ 1} = p_{1 2} - p_{2 1} \end{equation}\] so we ignore the components of the contingency table on the diagonal and the null hypothesis of no difference is that the off-diagonal components are equal. The hypothesis test with this null hypothesis and general alternative hypothesis is called McNemar’s test.

Let us try it on the data in Table 11.1 in Agresti.

mcnemar.test(y, correct = FALSE)
## 
##  McNemar's Chi-squared test
## 
## data:  y
## McNemar's chi-squared = 20.629, df = 1, p-value = 5.576e-06

This agrees with Section 11.1.3 in Agresti. He does not recommend R’s default correct = TRUE so we don’t either.

5 The Exact Competitor to McNemar’s Test

In McNemar’s test we are only interested in the off-diagonal components of the contingency table (whether they are equal). The theory of optimal tests taught in PhD level mathematical statistics and the basis of the fuzzy tests discussed by Geyer and Meeden says that, if do a conditional hypothesis test in which the test statistic is one of the off-diagonal components and we condition on the sum of the off-diagonal components, then the exact sampling distribution of the test statistic is \(\text{Bin}(y_{1 2} + y_{2 1}, 1 / 2)\) under the null hypothesis of no difference. Moreover, this test if one-tailed is uniformly most powerful (UMP) among all tests and if two-tailed is uniformly most powerful unbiased (UMPU) among all tests. So these are the best possible tests.

Although the tests are conditional, the UMP or UMPU property is unconditional. These are the best possible tests unconditionally.

Hence we know how to do hypothesis tests, even fuzzy tests, using the binomial distribution.

5.1 Conventional, Non-Fuzzy

These tests are conservative-exact rather than exact-exact as fuzzy tests are. We can not have exact significance level \(\alpha\) for any \(\alpha\) as explained in the notes about fuzzy tests. The test that rejects the null hypothesis whenever \(P \le \alpha\) is conservative. Its rejection probability is no greater than \(\alpha\).

Let’s try it on the voting data.

# lower-tailed
pbinom(y[1 , 2], y[1, 2] + y[2, 1], 1 / 2)
## [1] 2.926978e-06
# upper-tailed
pbinom(y[2 , 1], y[1, 2] + y[2, 1], 1 / 2)
## [1] 0.9999992
# two-tailed
2 * min(pbinom(y[1 , 2], y[1, 2] + y[2, 1], 1 / 2),
    pbinom(y[2 , 1], y[1, 2] + y[2, 1], 1 / 2))
## [1] 5.853956e-06

In the two-tailed case we get almost the same \(P\)-value at the asymptotic calculation done by R function mcnemar.test.

Moreover, \(P = 6\times 10^{-6}\) is so small that we aren’t actually worried about the test not being exact. But we keep going and do an exact fuzzy test.

5.2 Fuzzy

So these tests follow the notes on fuzzy hypothesis tests except for being about McNemar type data.

The lower-tailed test (which we know from the above is the one that gives a statistically significant result) is uniformly distributed on the interval

pbinom(c(y[1, 2] - 1, y[1 , 2]), y[1, 2] + y[2, 1], 1 / 2)
## [1] 8.262599e-07 2.926978e-06

And, because we have a symmetric distribution of the test statistic under the null hypothesis (success probability \(1 / 2\)), the \(P\)-value for two tails is twice the \(P\)-value for one-tail. So we are uniformly distributed on the interval

(2 * pbinom(c(y[1, 2] - 1, y[1 , 2]), y[1, 2] + y[2, 1], 1 / 2)) |> pmin(1)
## [1] 1.652520e-06 5.853956e-06

The funny pmin(1) at the end is to keep the range of the fuzzy \(P\)-value less than one in the corner case that the test statistic is exactly half of the conditioning count, that is, \(y_{1 2} = y_{2 1}\).

6 Confidence Interval Dual to McNemar’s Test

Agresti Section 11.1.1 suggests we use \(p_{1 2} - p_{2 1}\) to estimate \(\pi_{1 2} - \pi_{2 1}\). The variance is \[\begin{align*} \mathop{\rm var}(p_{1 2} - p_{2 1}) & = \frac{\mathop{\rm var}(y_{1 2} - y_{2 1})}{n^2} \\ & = \frac{\mathop{\rm var}(y_{1 2}) + \mathop{\rm var}(y_{2 1}) - 2 \mathop{\rm cov}(y_{1 2}, y_{2 1})}{n^2} \\ & = \frac{\pi_{1 2} (1 - \pi_{1 2}) + \pi_{2 1} (1 - \pi_{2 1}) + 2 \pi_{1 2} \pi_{2 1}}{n} \\ & = \frac{\pi_{1 2} + \pi_{2 1} - (\pi_{1 2} - \pi_{2 1})^2}{n} \end{align*}\] This agrees with Agresti’s formula.

pi.hat <- y / sum(y)
pi.hat
##             Election2008
## Election2004  Democrat Republican
##   Democrat   0.4041570  0.0369515
##   Republican 0.1247113  0.4341801
v <- pi.hat[1, 2] + pi.hat[2, 1] - (pi.hat[1, 2] - pi.hat[2, 1])^2
v
## [1] 0.153961
se <- sqrt(v / sum(y))
conf.level <- 0.95
crit <- qnorm((1 + conf.level) / 2)
crit
## [1] 1.959964
pi.hat[1, 2] - pi.hat[2, 1] + c(-1, 1) * crit * se
## [1] -0.12471791 -0.05080172

That is the Wald confidence interval for \(\pi_{1 2} - \pi_{2 1}\). Agresti agrees but has the opposite sign because he does the subtraction the opposite way (there is no right way, either is correct as long as one is not confused). Agresti also mentions a score interval, and a likelihood interval would also be possible. But since we have more exact conditional procedures, we leave this here.

7 The Exact Competitor to the Confidence Interval Dual to McNemar’s Test

The conditional distribution of \(y_{1 2}\) given \(y_{1 2} + y_{2 1}\) is \[ \text{Bin}\left( y_{1 2} + y_{2 1}, \frac{\pi_{1 2}}{\pi_{1 2} + \pi_{2 1}}\right) \] so the conservative-exact (conventional) or exact-exact (fuzzy) procedure must involve the parameter of this distribution. Call it \(\psi\).

7.1 Conventional, Non-Fuzzy

So we know how to make Wald, Wilks, or Rao intervals for this parameter based on this conditional distribution. For example, the Wald interval is

psi.hat <- pi.hat[1, 2] / (pi.hat[1, 2] + pi.hat[2, 1])
psi.hat
## [1] 0.2285714
se.psi.hat <- sqrt(psi.hat * (1 - psi.hat) / (y[1, 2] + y[2, 1]))
psi.hat + c(-1, 1) * crit * se.psi.hat
## [1] 0.1302025 0.3269403

7.2 Fuzzy

A fuzzy confidence interval is a graph showing the membership function.

fci.binom(y[1, 2], y[1, 2] + y[2, 1], alpha = 1 - conf.level, flat = 8)
## 95 percent fuzzy confidence interval
## core is [0.1465, 0.3279]
## support is (0.1349, 0.3438)

The interesting bit here is that not only do we get an exact fuzzy interval but we have to change the parameter from \(\pi_{1 2} - \pi_{2 1}\) to \(\pi_{1 2} / (\pi_{1 2} + \pi_{2 1})\) to make an exact confidence interval. This is the fuzzy confidence interval dual to the UMPU (randomized or fuzzy) competitor to McNemar’s test.

8 Generalizing Beyond Two-by-Two

We are now going to skip Section 11.2 in Agresti, going to Sections 11.3 and 11.4.

We are still working with square tables, just with more than two categories in each dimension. \(\pi_{a b}\) is still the probability that one individual is classified in category \(a\) before and in category \(b\) after (or left and right, or whatever the matched pairs are).

But now we have no equation like (\({*}\)) in the section about McNemar’s test above to simplify interpretation.

We are still interested in the marginal totals, which Agresti writes as \(\pi_{a +}\) for the before marginal totals and \(\pi_{+ b}\) for the after marginal totals.

But just constraining the marginal totals to be equal does not give a log-linear model, hence not a regular full exponential family. Thus we look at some other models first.

8.1 Data

Table 11.5 in Agresti has data on migration.

data(table_11.5)
foo <- xtabs(count ~ row + column, data = table_11.5)
names(dimnames(foo)) <- c("residence age 16", "residence 2010")
foo
##                 residence 2010
## residence age 16  mw  ne   s   w
##               mw 414  10  50  40
##               ne  15 266  61  28
##               s   22   8 578  22
##               w    6   7  27 301

8.2 Symmetry

If we write the saturated model with canonical statistic and canonical parameter as matrices with canonical bilinear form \[ \langle y, \theta \rangle = \sum\nolimits_a \sum\nolimits_b y_{a b} \theta_{a b} \] and then impose symmetry \[ \theta_{a b} = \theta_{b a}, \qquad \text{for all $a$ and $b$}, \] then we can rewrite the canonical bilinear form as \[ \langle y, \theta \rangle = \left( \sum_{a} y_{a a} \theta_{a a} \right) + \left( \sum_{\substack{(a, b) \\ a < b}} (y_{a b} + y_{b a}) \theta_{a b} \right) \] Symmetry of the canonical parameter implies symmetry for the mean value parameter because one is the log of the other times a constant. Hence the likelihood equations are (by observed equals expected) \[ \hat{\pi}_{a b} = \frac{y_{a b} + y_{b a}}{2 n} \] where \(n\) is the cell total.

pi.hat.sym <- (foo + t(foo)) / (2 * sum(foo))
pi.hat.sym
##                 residence 2010
## residence age 16          mw          ne           s           w
##               mw 0.223180593 0.006738544 0.019407008 0.012398922
##               ne 0.006738544 0.143396226 0.018598383 0.009433962
##               s  0.019407008 0.018598383 0.311590296 0.013207547
##               w  0.012398922 0.009433962 0.013207547 0.162264151

The likelihood ratio test statistic with symmetry as the null hypothesis and the general model as the alternative hypothesis is

lrt <- sum(foo * log(foo / sum(foo) / pi.hat.sym))
lrt
## [1] 50.24137

If we assume Poisson sampling the degrees of freedom is \(d^2\) for the general model (alternative hypothesis) and \(d (d + 1) / 2\) for the symmetry model (null hypothesis), where \(d\) is the dimension, so the \(P\)-value is

d <- nrow(foo)
df <- d^2 - d * (d + 1) / 2
pchisq(lrt, df, lower.tail = FALSE)
## [1] 4.205337e-09

Here we are thinking of this test as a goodness of fit test. We want symmetry to fit the data. It clearly does not. This model is way too restrictive. We need some more interesting models.

8.3 Quasi-Symmetry

If we write our model as a general hierarchical model with \[\begin{equation} \tag{${*}{*}$} \theta_{a b} = \alpha + \beta_a + \gamma_b + \delta_{a b} \end{equation}\] and then impose symmetry only on the interaction term but not on the main effects terms (that is, we do not impose \(\beta_a = \gamma_a\) for all \(a\)), then the resulting model is called quasi-symmetry (Agresti, Section 11.4.2).

AFAIK there is no way to make the R formula language describe this model so R function glm can fit it. But we can make a model matrix.

But in order to get an identifiable model, we must first consider identifiability because we cannot have R function glm figure that out for us.

In (\({*}{*}\))

  • \(\delta_{+ b}\) is confused with \(\gamma_b\),

  • \(\delta_{a +}\) is confused with \(\beta_a\),

  • \(\gamma_{+}\) is confused with \(\alpha\), and

  • \(\beta_{+}\) is also confused with \(\alpha\).

Hence we can gain identifiability by setting one row and one column of \(\delta_{a b}\) to zero and setting one component of each of \(\beta_a\) and \(\gamma_b\) to zero. And setting them equal to zero is equivalent to leaving the corresponding dummy variables out of the model matrix. We will leave out the terms for mw.

names(table_11.5)
## [1] "row"    "column" "count"
resp <- table_11.5$count
modmat <- rep(1, length(resp))
modmat <- cbind(modmat)
colnames(modmat) <- "intercept"
for (i in setdiff(table_11.5$row, "mw")) {
    dumbo <- as.integer(table_11.5$row == i)
    dumbo <- cbind(dumbo)
    colnames(dumbo) <- paste0("row", i)
    modmat <- cbind(modmat, dumbo)
}
for (i in setdiff(table_11.5$col, "mw")) {
    dumbo <- as.integer(table_11.5$col == i)
    dumbo <- cbind(dumbo)
    colnames(dumbo) <- paste0("col", i)
    modmat <- cbind(modmat, dumbo)
}
for (i in setdiff(table_11.5$row, "mw"))
    for (j in setdiff(table_11.5$col, "mw"))
        if (i <= j) {
            dumbo <- as.integer(table_11.5$row == i & table_11.5$col == j |
                table_11.5$row == j & table_11.5$col == i)
            dumbo <- cbind(dumbo)
            colnames(dumbo) <- paste0(i, ":", j)
            modmat <- cbind(modmat, dumbo)
        }
modmat
##       intercept rowne rows roww colne cols colw ne:ne ne:s ne:w s:s s:w w:w
##  [1,]         1     1    0    0     1    0    0     1    0    0   0   0   0
##  [2,]         1     1    0    0     0    0    0     0    0    0   0   0   0
##  [3,]         1     1    0    0     0    1    0     0    1    0   0   0   0
##  [4,]         1     1    0    0     0    0    1     0    0    1   0   0   0
##  [5,]         1     0    0    0     1    0    0     0    0    0   0   0   0
##  [6,]         1     0    0    0     0    0    0     0    0    0   0   0   0
##  [7,]         1     0    0    0     0    1    0     0    0    0   0   0   0
##  [8,]         1     0    0    0     0    0    1     0    0    0   0   0   0
##  [9,]         1     0    1    0     1    0    0     0    1    0   0   0   0
## [10,]         1     0    1    0     0    0    0     0    0    0   0   0   0
## [11,]         1     0    1    0     0    1    0     0    0    0   1   0   0
## [12,]         1     0    1    0     0    0    1     0    0    0   0   1   0
## [13,]         1     0    0    1     1    0    0     0    0    1   0   0   0
## [14,]         1     0    0    1     0    0    0     0    0    0   0   0   0
## [15,]         1     0    0    1     0    1    0     0    0    0   0   1   0
## [16,]         1     0    0    1     0    0    1     0    0    0   0   0   1

Inelegant. But works.

moo <- glm.fit(modmat, resp, family = poisson())$fitted
moo <- xtabs(moo ~ row + column, data = table_11.5)
names(dimnames(moo)) <- c("residence age 16", "residence 2010")
moo
##                 residence 2010
## residence age 16         mw         ne          s          w
##               mw 414.000000   9.286764  55.245544  35.467692
##               ne  15.713236 266.000000  58.512325  29.774439
##               s   16.754456  10.487675 578.000000  24.757869
##               w   10.532308   5.225561  24.242131 301.000000
foo
##                 residence 2010
## residence age 16  mw  ne   s   w
##               mw 414  10  50  40
##               ne  15 266  61  28
##               s   22   8 578  22
##               w    6   7  27 301

According to equation (11.20) in Agresti, the observed and MLE expected cell counts should satisfy the following.

all.equal(rowSums(moo), rowSums(foo))
## [1] TRUE
all.equal(colSums(moo), colSums(foo))
## [1] TRUE
all.equal(moo + t(moo), foo + t(foo), check.attributes = FALSE)
## [1] TRUE

And now for the big test: how well does this model fit the data. We do LRT again.

lrt <- sum(foo * log(foo / moo))
df <- d^2 - ncol(modmat)
pchisq(lrt, df, lower.tail = FALSE)
## [1] 0.3218021

So this model fits the data well. The general hypothesis (alternative) fits no better than can be expected for having more degrees of freedom.

And this is a log-linear model, hence a regular full exponential family with all the nice properties of such.

So now that we have found this model fits, it is interesting? And for that we apply our philosophy of exponential family models. They make sense if the canonical statistics make sense. So what are they? \[\begin{align*} \langle y, \theta \rangle & = \sum\nolimits_{a b} y_{a b} \theta_{a b} \\ & = \sum\nolimits_{a b} y_{a b} (\alpha + \beta_a + \gamma_b + \delta_{a b}) \\ & = y_{+ +} \alpha + \left(\sum\nolimits_a y_{a +} \beta_a \right) + \left(\sum\nolimits_a y_{+ b} \gamma_b \right) + \left(\sum\nolimits_a y_{a a} \delta_{a a} \right) + \left(\sum_{\substack{a b \\ a < b}} (y_{a b} + y_{b a}) \delta_{a b} \right) \end{align*}\] Thus the canonical sufficient statistics are

  • \(y_{+ +}\) which is needed only if we have Poisson sampling (if multinomial sampling, it is the known sample size, hence not random),

  • \(y_{a +}\) the row margins,

  • \(y_{+ b}\) the column margins,

  • \((y_{a b} + y_{b a}) / 2\), the symmetrization of the response (also called the symmetric part).

And we also see that the model fitting equations (Agresti, equation (11.20) are just observed equals expected for the quasi-symmetry canonical statistics.

Margins make scientific sense as well as statistical sense. But does symmetrization? That is a scientific question, not a statistical one. So we do not have an answer for that. We leave it to the scientists in each particular application.

We can, however, apply the principles for explaining exponential families. We do have both marginals as part of the canonical sufficient statistic vector. Thus the model does adjust (observed equals expected) for the populations of the regions being different (proportionally) at the two times. Then we also have the symmetrization of the data matrix as part of the canonical sufficient statistic. This means the model does not keep track of moves from region A to region B or vice versa. It only keeps track of moves between each pair of regions, not which way the moves go. This does not mean that the model assumes moves in both directions are equally probable. The different margins assure that does not happen. To apply our principles, sufficiency says which way the moves go does not matter scientifically, and maximum entropy says which way the moves go is a random as it is possible to be (in the sense of maximum entropy) subject to observed equals expected for the canonical statistics.

Note also that the symmetry model is a submodel of the quasi-symmetry model obtained by setting \(\beta_a = \gamma_a\) for all \(a\).

8.4 Marginal Homogeneity

So for our next trick, we drop symmetry (quasi- or otherwise) altogether. We just require the margins agree.

Unfortunately, this is not a log-linear model, so R function glm is useless. We will have to use the general principles of maximum likelihood estimation (code up the log likelihood and maximize it). But now we have to maximize it subject to constraints (equal margins). As we did before (calculating likelihood-based confidence intervals we will use R function auglag in CRAN package alabama for constrained optimization.

Here, since the constraints are linear in the mean-value parameters (or the usual, cell probability, parameters) we will use those parameters. The multinomial log-likelihood is \[ l(\pi) = \sum\nolimits_{a b} y_{a b} \log(\pi_{a b}) \] and the constraints are that the components of \(\pi_{a b}\) are nonnegative and sum to one and (for the hypothesis of marginal homogeneity) \(\pi_{a +} = \pi_{+ a}\) for all \(a\).

So let’s try it. We need three functions: minus log likelihood, inequality constraint function, and equality constraint function. To avoid global variables, we use a function factory. To properly communicate between the three functions we use a function factory to make their function factories.

combo.factory <- function(y) {
    stopifnot(is.numeric(y))
    stopifnot(is.finite(y))
    stopifnot(is.matrix(y))
    stopifnot(nrow(y) == ncol(y))
    stopifnot(y == round(y))
    stopifnot(y >= 0)
    function(type = c("mlogl", "inequality", "equality")) {
        ty <- match.arg(type)
        if (ty == "mlogl")
            return(function(pi) {
                stopifnot(is.numeric(pi))
                stopifnot(is.finite(pi))
                stopifnot(length(pi) == length(y))
                pi <- matrix(pi, nrow(y), ncol(y))
                # other constraints left to constraint function
                # do not return NaN for log(pi) when pi < 0
                foo <- suppressWarnings(log(pi))
                foo[is.na(foo)] <- -Inf
                sum(- y * foo)
            })
        if (ty == "inequality")
            return(function(pi) {
                stopifnot(is.numeric(pi))
                stopifnot(is.finite(pi))
                stopifnot(length(pi) == length(y))
                pi <- matrix(pi, nrow(y), ncol(y))
                return(pi)
            })
        if (ty == "equality")
            return(function(pi) {
                stopifnot(is.numeric(pi))
                stopifnot(is.finite(pi))
                stopifnot(length(pi) == length(y))
                pi <- matrix(pi, nrow(y), ncol(y))
                foo <- rowSums(pi) - colSums(pi)
                foo <- c(foo, sum(pi) - 1)
                return(foo)
            })
    }
}

The help for R function auglag says we do not need to start with feasible point, but might as well.

pi.start <- matrix(1 / length(foo), nrow(foo), ncol(foo))
pi.start
##        [,1]   [,2]   [,3]   [,4]
## [1,] 0.0625 0.0625 0.0625 0.0625
## [2,] 0.0625 0.0625 0.0625 0.0625
## [3,] 0.0625 0.0625 0.0625 0.0625
## [4,] 0.0625 0.0625 0.0625 0.0625
cout <- combo.factory(foo)
aout <- auglag(pi.start, fn = cout("mlogl"), hin = cout("in"),
    heq = cout("eq"), control.outer = list(trace = FALSE))
aout$convergence == 0
## [1] TRUE

Guess we have the solution. Anyway, it says it has the solution.

pi.hat <- aout$par
dimnames(pi.hat) <- dimnames(foo)
pi.hat
##                 residence 2010
## residence age 16          mw          ne          s           w
##               mw 0.223196268 0.006621942 0.01769656 0.014119457
##               ne 0.006793672 0.143521385 0.01919783 0.008787421
##               s  0.024830708 0.014858048 0.31150302 0.011798066
##               w  0.006813569 0.013298933 0.01459244 0.162370682
rowSums(pi.hat) - colSums(pi.hat)
##            mw            ne             s             w 
##  5.863252e-09  3.033819e-09 -3.925868e-09 -4.971203e-09

Constraints pretty close to zero.

Table 11.5 in Agresti also appears to have model fits.

mu.hat <- sum(foo) * pi.hat
round(mu.hat, 1)
##                 residence 2010
## residence age 16    mw    ne     s     w
##               mw 414.0  12.3  32.8  26.2
##               ne  12.6 266.2  35.6  16.3
##               s   46.1  27.6 577.8  21.9
##               w   12.6  24.7  27.1 301.2

These fitted values seem to agree with those in Agresti.

So finally for the likelihood ratio tests. This is a supermodel of symmetry and a submodel of general. Versus general first.

lrt <- sum(foo * log(foo / sum(foo) / pi.hat))
pchisq(lrt, df = d - 1, lower.tail = FALSE)
## [1] 3.798478e-10

Marginal homogeneity does not fit these data. Clearly, people are moving at different rates into and out of some region or regions. But, of course, this test says nothing about which region or regions.

9 Independence and Quasi-Independence

In Section 11.4.4 Agresti introduces another model: quasi-independence. We will cover that but first we do independence from intro stats.

chisq.test(foo)
## 
##  Pearson's Chi-squared test
## 
## data:  foo
## X-squared = 3433.4, df = 9, p-value < 2.2e-16

That does not fit. Not close to complete random movement (where they were in 2010 completely unrelated to where they were at age 16).

The quasi-independence is like the independence model except for extra mass on the diagonal. The model equation is \[\begin{equation} \tag{${*}{*}{*}$} \theta_{a b} = \alpha + \beta_a + \gamma_b + \delta_a I(a = b) \end{equation}\] (Agresti, equation (11.26)) where \(I(\,\cdot\,)\) is a function that maps false to zero and true to one. We need to turn that into what R calls a factor variable.

diag <- table_11.5$row
diag[table_11.5$row != table_11.5$column] <- "off"
diag
##  [1] "ne"  "off" "off" "off" "off" "mw"  "off" "off" "off" "off" "s"   "off"
## [13] "off" "off" "off" "w"
qiout <- glm(count ~ row + column + diag, data = table_11.5, family = poisson)
moo.qi <- qiout$fitted
moo.qi <- xtabs(moo.qi ~ row + column, data = table_11.5)
names(dimnames(moo.qi)) <- names(dimnames(foo))
moo.qi
##                 residence 2010
## residence age 16         mw         ne          s          w
##               mw 414.000000  10.876876  56.976620  32.146505
##               ne  18.295896 266.000000  54.790832  30.913272
##               s   15.944463   9.115313 578.000000  26.940223
##               w    8.759640   5.007811  26.232548 301.000000

And another LRT against the general hypothesis.

lrt <- sum(foo * log(foo / moo.qi))
df <- length(foo) - length(qiout$coefficients)
pchisq(lrt, df, lower.tail = FALSE)
## [1] 0.4661752

So we now have two models that fit the data: quasi-symmetry and quasi-independence. Can we say anything about which fits better? Comparing our equations \(({*}{*})\) in the section about quasi-symmetry and \(({*}{*}{*})\) in this section, we see that the quasi-independence model is a special case of the quasi-symmetry model. The former just constrains all of the \(\delta_{a b}\) for \(a \neq b\) to be zero.

So let’s do this test.

lrt <- sum(foo * log(moo / moo.qi))
df <- ncol(modmat) - length(qiout$coefficients)
lrt
## [1] 1.111362
df
## [1] 2

Agresti says the quasi-independence model has \(d\) more parameters than the independence model, where \(d\) is the dimension. And we know from the chapter zero notes that the independence model has \(2 d - 1\) identifiable parameters, so that gives quasi-independence \(3 d - 1\) identifiable parameters. Check if we are right.

length(qiout$coefficients)
## [1] 11
3 * nrow(foo) - 1
## [1] 11

So any worries about the above model fit seem misplaced. The R formula mini-language and R function glm seem to have nailed it. Hence our \(P\)-value calculation

pchisq(lrt, df = df, lower.tail = FALSE)
## [1] 0.5736815

says we clearly have no reason to reject the null hypothesis of quasi-independence. So we cannot tell whether quasi-symmetry or quasi-independence fits these data better.

If we want to pick the smaller model (quasi-independence) on the basis of parsimony, we can do that. But statistics isn’t telling us to do that.

We also have a clear scientific interpretation of quasi-independence. Ignoring the non-movers (the cells on the diagonal, which both models fit perfectly, the counts in the cells on the diagonal being sufficient statistics for both, so observed equals expected for them), quasi-independence says that where people move (inside the United States) is independent of where the were at age 16.

So quasi-independence has a clearer scientific interpretation than quasi-symmetry AFAICS.

10 Higher Dimensional Tables

Agresti (still in Chapter 11) discusses analogs of independence, symmetry, quasi-symmetry, and quasi-independence for three-way, four-way, and so forth, tables (so for data on matched triplets, quadruplets, etc.).

The basic idea is that we may impose any of these (independence, symmetry, quasi-symmetry, or quasi-independence) on some margins of a table but not on other margins, especially on lower-dimensional margins, but not higher-dimensional margins. But Agresti says

This hierarchy is mathematically attractive. However, the higher-order symmetries are usually too restrictive to fit well in practice.

So he does no examples and we won’t either.