Warp-breaks data

Section of Marginal models, and two responses at once, part of Correlated data

Warp-breaks data

Type	Section
Taught in	Marginal models, and two responses at once
Part	Correlated data

Yarn was spun at two factories, at low, medium, and high tension, and the number of breaks in a fixed length was recorded, nine replicates of each tension at each factory. Breaks are a count. Yarn from the same factory shares machines, material, and a standard, so the factory is the cluster. Two clusters is a small number for a sandwich. The model is still fitted. The robust column is the one the theory prefers, and it is not a precise standard error from a large number of factories.

An exchangeable working correlation says every pair of looms from the same factory has the same correlation. With a Poisson mean and a log link,

need <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, repos = "https://cloud.r-project.org")
  }
  suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}
need("gee")
data(warpbreaks)
warpbreaks$factory <- warpbreaks$wool
warpbreaks <- warpbreaks[order(warpbreaks$factory), ]
exchangeable <- gee(breaks ~ tension, id = factory, data = warpbreaks,
                    family = poisson, corstr = "exchangeable")
summary(exchangeable)
ar1 <- gee(breaks ~ tension, id = factory, data = warpbreaks,
           family = poisson, corstr = "AR-M", Mv = 1)
summary(ar1)
ar1$working.correlation[1:4, 1:4]
gaussian_ar <- gee(breaks ~ tension, id = factory, data = warpbreaks,
                   family = gaussian(link = "log"), corstr = "AR-M", Mv = 1)
summary(gaussian_ar)

            Estimate   Naive SE   Robust SE
(Intercept)   3.594     0.091       0.159
tensionM     -0.321     0.128       0.223
tensionH     -0.518     0.136       0.064

The scale estimate is 4.60, far from 1: these counts are more variable than a Poisson. The working correlation is about 0.021. Low tension, the baseline, has fitted mean \exp(3.594)\approx 36 breaks. Medium tension multiplies that mean by \exp(-0.321)=0.73. High tension multiplies it by \exp(-0.518)=0.60. On the naive standard errors, medium tension is detected (z=-2.51). On the robust standard errors it is not (z=-1.44). High tension is detected either way, and the robust standard error is about half the naive one. The two columns disagree, which is evidence that the Poisson variance and the exchangeable correlation are not the truth. The robust column is the one that remains consistent if the mean is right.

The looms were recorded in the order they came off the machine. Then exchangeable correlation is the wrong working model. An AR(1) working correlation sets \mathrm{Corr}(Y_t,Y_{t+k})=\gamma^{|k|}. The lag-1 value is 0.041, then 0.041^2 at lag 2 and 0.041^3 at lag 3. The coefficients barely move: intercept 3.595, medium -0.325, high -0.518. The naive standard errors move, because the working model moved. The robust standard errors stay where they were, to two decimals. That is the theorem, visible in one table. A third fit keeps the log mean and the AR(1) working correlation, and replaces the Poisson variance by a constant variance, a Gaussian family with a log link. A constant-variance model for counts is a poor variance story. The coefficients are again 3.595, -0.324, -0.518, and the robust standard errors again match the previous robust column. The scale estimate, 141, is a variance on the count scale and is not a quantity to interpret as a Poisson dispersion.

Key idea

If the mean is the right function of the covariates, the estimates of \boldsymbol{\beta} and the sandwich standard errors settle in the same place whatever working correlation the fit starts from. The naive standard errors move. A large gap between the two columns means the working variance was wrong. It does not, by itself, mean the robust column is wrong.

The gap here is large for every working model fitted above, so none of them is close to the true covariance. The robust standard errors are still the ones the asymptotic result licenses. They are not a substitute for having more than two factories.

Warp-breaks data

Used in

See also