How to best find standard error *across* linear regression fits?

Question

So I have a scenario where there are $n = 8$ subjects, which are observed at 20 time points and having heteroscedasticity in their response. For example, consider the following:

num_datasets = 8;

x = [1:20]';

%define matrix for the response for 8 different datasets
Y = repmat(x,1,8) * nan;

for i = 1:size(X,2)
    Y(:,i) = 2*x + unifrnd(3,8)*randn(size(x));
end

So clearly each observation/subject has the same linear model relating their response ($y$) to the regressor ($x$), but the amounts/sources of noise vary between subject. Now, I know that the standard error for the linear regression fit has the form:

$$\sigma\sqrt{\frac{1}{n}+ \frac{(x^*-\bar x)^2}{\sum_{i=1}^n (x_i-\bar{x})^2} }$$

where $\sigma$ represents the standard deviation of the residuals of the fit, $n$ represents the number of samples in the observation (in my example above this would be 20, not 8), $(x^* - \bar x)$ represents the distance of each $x_i$ sample from the mean (which is why the standard error increases hyperbolically as you deviate from the mean), and then ${\sum_{i=1}^n (x_i-\bar{x})^2}$ is simply the variance in $x$.

However, if I interpret this equation correctly, I think this gives the standard error across the dimension of $x$, and doesn't directly tell me the standard error across subjects. In other words, I suspect it wouldn't be a good idea to use this formula for each subject and then take the mean standard error (please correct me if I am wrong). So I have 2 questions:

1. What would be the best way to calculate the standard error across subjects? Would it simply be to perform the fit for each subject, and then take the standard deviation of the fits?

2. What would the shape of the standard error of the fit look like, and what is the intuition behind that? Would it still be hyperbolic? I don't think it would, but actually really not sure.

Answer 1

It might be best to think about this situation in terms of meta-analysis: putting together information from several studies to estimate the model for the underlying population. Studies are combined by weighting them according to information they provide, typically inversely weighting each study by the variance of its estimates.

You can think about your case as representing 8 different "studies" (8 different subjects), with each having values of $y$ measured at 20 values of $x$. We assume that within each subject the standard assumptions of linear regression hold, in particular that observations are uncorrelated and the variance of $y$ about the regression is independent of the value of $x$. Unlike many practical meta-analyses that depend on reported summaries of results for each of several studies, you still have the individual data for each "study."

So if you want a model for the underlying population, one simple way to proceed would be to do each of the 8 individual regressions and determine the value of estimated residual variance $\hat\sigma_j^2$ for each subject $j$.* Then re-weight each individual data point inversely to that estimated variance for the corresponding subject, and perform a weighted least-squares regression over all the 160 data points.

What you call the "hyperbolic" shape of error in $\hat y$ for new predictions as a function of $x$ will be the same. It comes from the uncertainty in the estimate of the slope in the regression. The error is smallest ($\sigma/\sqrt{n}$) at the mean value of $x$, then increases with that shape simply because you are uncertain how quickly $y$ changes with $x$ as you move farther away from $\bar x$. Linear regression with uncorrelated observations weighted by their variances provides the best linear unbiased estimators (BLUE) of the regression coefficient,** and with the weighted regression combining all cases you now have an $n$ value of 160. So the width of that uncertainty area will tend to be minimized.

Some caution is in order, as the NIST page puts it:

The biggest disadvantage of weighted least squares, which many people are not aware of, is probably the fact that the theory behind this method is based on the assumption that the weights are known exactly. This is almost never the case in real applications, of course, so estimated weights must be used instead.

And as @cardinal put it:

Learning a variance is hard.

For a normal distribution with variance $\sigma^2$, the variance of a variance estimate $\hat\sigma^2$ from $n$ observations is $2\sigma^4/(n-1)$. So unless you have many data points and a reason to believe that there are substantial differences in the true $\sigma_j^2$ values among the subjects $j$ there might not be much benefit to this weighting approach.

The above assumes that all subjects have the same slopes and intercepts for the relationship between $y$ and $x$. One might interpret your suggestion to "take the standard deviation of the fits" as meaning that you expect true differences among subjects in these parameter values. In that case you could obtain estimates for the variances of intercepts and slopes among subjects with a mixed model. Weighting of the individual data points could still be done.

---

*I believe that there is a way to estimate all the within-subject variances and the shared regression coefficients in a single model, but I don't recall immediately what that is. It would probably require an iterative or maximum-likelihood approach. This is simple approach gets to the substance of your question.

**See the page linked on weighted regression. This assumes that the variances are known.

Answer 2

Question 1. What would be the best way to calculate the standard error across subjects? Would it simply be to perform the fit for each subject, and then take the standard deviation of the fits?

Option 1: Use weighted least squares. The Gauss Markov theorem tells us that the inverse variance weighted standard error estimator will be the Best Linear Unbiased Estimator (BLUE). Note that although the mean model is correct, and consequently the unweighted estimate is unbiased, there is added efficiency of using the iterative generalized least squares estimator to provide a better estimate of the residuals. It helps to identify the appropriate degrees of freedom for the intracluster variance. For reference, I've included the two-stage estimate, but I have trouble identifying the correct degree of freedom correction.

One interesting result that I'm working on is the idea that off-the-box software with options for intracluster correlation can provide consistent estimates of heteroscedasticity. That is, regardless of whether a sample is highly intracorrelated or highly variable, the net effect is a downweighting of that sample, and so the same optimal standard error can be obtained in either case.

Using unweighted residuals to estimate cluster variance, I'm finding is that it's hard to identify the appropriate degree of freedom for the intracluster variance estimate. I'm adding my code below for others to verify. $n-1$ is too conservative, and $n-2$ is too conservative.

Option 2: Use the sandwich variance estimator (heteroscedasticity consistent) or bootstrap.

Question 2: What would the shape of the standard error of the fit look like, and what is the intuition behind that? Would it still be hyperbolic? I don't think it would, but actually really not sure.

The limiting distribution of the error distribution is still normal provided that the sample "grows faster" in terms of number of time points rather than in the number of subjects, or at least such that subject level heteroscedasticity is somewhat bounded. The intuition is that it is a result of the Lyapunov central limit theorem.

require(gee)
`%covers%` <- function(x, y) x[1] < y & y < x[2]
sse.df <-function(x, df=1) {
  sum({x-mean(x)}^2)/{length(x)-df}
}
confint.gee <- function (object, parm, level = 0.95, ...) 
{
  cf <- coef(object)
  pnames <- names(cf)
  if (missing(parm)) 
    parm <- pnames
  else if (is.numeric(parm)) 
    parm <- pnames[parm]
  a <- (1 - level)/2
  a <- c(a, 1 - a)
  # pct <- format.perc(a, 3)
  pct <- paste0(formatC(100*a, format='f', digits=1), '%')
  fac <- qnorm(a)
  ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, 
                                                             pct))
  # ses <- sqrt(diag(vcov(object)))[parm]
  ses <- sqrt(diag(object$robust.variance))[parm]
  ci[] <- cf[parm] + ses %o% fac
  ci
}


do.one <- function() {
  s1 <- 1
  s2 <- 1
  nc <- 8
  nt <- 20
  i <- rep(1:8, each=nt)
  
  e <- rnorm(nc, 0, s1)[i] + rnorm(nc*nt, 0, s2)
  
  x <- rep(seq(-3, 3, length.out = nt), times=nc)
  y <- 2*x + e
  
  r <- lm.fit(cbind(1,x), y)$residuals
  
  
  wls <- lm(y ~ x, weights=rep(1/tapply(r^2, i, sse.df, df=1), each=nt))
  gls <- gls(y ~ x, correlation=corCompSymm(form=~1|i))
  
  gee <- gee(y ~ x, id = i)
  
  c( ## coverage of 80% CIs
    confint(wls, parm='x', level = .8) %covers% 2,
    confint(gee, parm='x', level = .8) %covers% 2,
    confint(gls, parm='x', level= 0.8) %covers% 2,
    vcov(wls)[2,2]^.5,
    gee$robust.variance[2,2]^.5,
    vcov(gls)[2,2]^.5
  )
}

set.seed(123)
out <- replicate(500, do.one())

## 80% coverage of CIs
rowMeans(out[1:3, ])

par(mfrow=c(1,3))
hist(out[4, ], xlab='Sigma two-pass', main='')
hist(out[5, ], xlab='Sigma GEE', main='')
hist(out[6, ], xlab='Sigma GLS', main='')

Gives us coverage of 70% for 2 degree of freedom WLS and 74% for the GEE. and 82.54% for GLS. Histograms of the standard error estimates show a closely normal distribution in all cases.