10
$\begingroup$

For sake of a simple example assume that there are two linear regression models

  • Model 1 has three predictors, x1a, x2b, and x2c
  • Model 2 has three predictors from model 1 and two additional predictors x2a and x2b

There is a population regression equation where the population variance explained is $\rho^2_{(1)}$ for Model 1 and $\rho^2_{(2)}$ for Model 2. The incremental variance explained by Model 2 in the population is $\Delta\rho^2 = \rho^2_{(2)} - \rho^2_{(1)}$

I am interested in obtaining standard errors and confidence intervals for an estimator of $\Delta\rho^2$. While the example involves 3 and 2 predictors respectively, my research interest concerns a wide range of different numbers of predictors (e.g., 5 and 30). My first thought was to use $\Delta r^2_{adj} = r^2_{adj(2)} - r^2_{adj(1)}$ as an estimator and bootstrap it, but I wasn't sure whether this would be appropriate.

Questions

  • Is $\Delta r^2_{adj}$ a reasonable estimator of $\Delta \rho^2$?
  • How can a confidence interval be obtained for the population r-square change (i.e., $\Delta\rho^2$)?
  • Would bootstrapping $\Delta\rho^2$ be appropriate for confidence interval calculation?

Any references to simulations or the published literature would also be most welcome.

Example code

If it helps, I created a little simulation dataset in R which could be used to demonstrate an answer:

n <- 100
x <- data.frame(matrix(rnorm(n *5), ncol=5))
names(x) <- c('x1a', 'x1b', 'x1c', 'x2a', 'x2b')
beta <- c(1,2,3,1,2)
model2_rho_square <- .7
error_rho_square <- 1 - model2_rho_square
error_sd <- sqrt(error_rho_square / model2_rho_square* sum(beta^2))
model1_rho_square <- sum(beta[1:3]^2) / (sum(beta^2) + error_sd^2)
delta_rho_square <- model2_rho_square - model1_rho_square

x$y <- rnorm(n, beta[1] * x$x1a + beta[2] * x$x1b + beta[3] * x$x1c +
               beta[4] * x$x2a + beta[5] * x$x2b, error_sd)

c(delta_rho_square, model1_rho_square, model2_rho_square)
summary(lm(y~., data=x))$adj.r.square - 
        summary(lm(y~x1a + x1b + x1c, data=x))$adj.r.square

Reason for concern with bootstrap

I ran a bootstrap on some data with around 300 cases, and 5 predictors in the simple model and 30 predictors in the full model. While the sample estimate using adjusted r-square difference was 0.116, the boostrapped confidence interval were mostly larger CI95%(0.095 to 0.214) and the mean of the bootstraps was nowhere near the sample estimate. Rather the mean of the boostrapped samples appeared to be centred on the sample estimate of the difference between r-squares in the sample. This is despite the fact that I was using the sample adjusted r-squares to estimate the difference.

Interestingly, I tried an alternative way of computing $\Delta\rho^2$ as

  1. calculate sample r-square change
  2. adjust the sample r-square change using the standard adjusted r-square formula

When applied to the sample data this reduced the estimate of $\Delta \rho^2$ to .082 but the confidence intervals seemed appropriate for the the method I mentioned first, CI95% (.062, .179) with mean of .118.

Broadly, I'm concerned that bootstrapping assumes that the sample is the population, and therefore estimates that reduce for overfitting may not perform appropriately.

$\endgroup$
  • $\begingroup$ "However, I am concerned that bootstrapping such a population adjusted value might be problematic." -- why? $\endgroup$ – January Jun 28 '13 at 4:22
  • $\begingroup$ @January I've edited the question and tried to articulate my concern about bootstrapping with adjusted r-square. $\endgroup$ – Jeromy Anglim Jun 28 '13 at 5:54
  • $\begingroup$ What is the population R squared ? I've took a look at the definition given here but for me the variance $\sigma^2_y$ makes no sense because the $y_i$ are not identically distributed. $\endgroup$ – Stéphane Laurent Jun 28 '13 at 7:00
  • $\begingroup$ @StéphaneLaurent it is the percentage of variance explained in the population by the population regression equation. Or you could define it asymptotically as the proportion of variance explained in your sample as your sample size approaches infinity. See also this answer regarding unbiased estimates of population r-square. It is particularly relevant in psychology where we are often more interested in the true relationship rather than actually applying our estimated prediction equation. $\endgroup$ – Jeromy Anglim Jun 28 '13 at 8:16
  • 3
    $\begingroup$ An F-test can be thought of as test of the hypothesis $\Delta\rho^2 = 0$. Can that be used to derive the standard error and confidence interval you are looking for? $\endgroup$ – Maarten Buis Jun 28 '13 at 9:09
3
$\begingroup$

Population $R^2$

I'm firstly trying to understand the definition of the population R-squared.

Quoting your comment:

Or you could define it asymptotically as the proportion of variance explained in your sample as your sample size approaches infinity.

I think you mean this is the limit of the sample $R^2$ when one replicates the model infinitely many times (with the same predictors at each replicate).

So what is the formula for the asymptotic value of the sample $R^²$ ? Write your linear model $\boxed{Y=\mu+\sigma G}$ as in https://stats.stackexchange.com/a/58133/8402, and use the same notations as this link.
Then one can check that the sample $R^2$ goes to $\boxed{popR^2:=\dfrac{\lambda}{n+\lambda}}$ when one replicates the model $Y=\mu+\sigma G$ infinitely many times.

As example:

> ## design of the simple regression model lm(y~x0)
> n0 <- 10
> sigma <- 1
> x0 <- rnorm(n0, 1:n0, sigma)
> a <- 1; b <- 2 # intercept and slope
> params <- c(a,b)
> X <- model.matrix(~x0)
> Mu <- (X%*%params)[,1]
> 
> ## replicate this experiment k times 
> k <- 200
> y <- rep(Mu,k) + rnorm(k*n0)
> # the R-squared is:
> summary(lm(y~rep(x0,k)))$r.squared 
[1] 0.971057
> 
> # theoretical asymptotic R-squared:
> lambda0 <- crossprod(Mu-mean(Mu))/sigma^2
> lambda0/(lambda0+n0)
          [,1]
[1,] 0.9722689
> 
> # other approximation of the asymptotic R-squared for simple linear regression:
> 1-sigma^2/var(y)
[1] 0.9721834

Population $R^2$ of a submodel

Now assume the model is $\boxed{Y=\mu+\sigma G}$ with $H_1\colon\mu \in W_1$ and consider the submodel $H_0\colon \mu \in W_0$.

Then I said above that the population $R^2$ of model $H_1$ is $\boxed{popR^2_1:=\dfrac{\lambda_1}{n+\lambda_1}}$ where $\boxed{\lambda_1=\frac{{\Vert P_{Z_1} \mu\Vert}^2}{\sigma^2}}$ and $Z_1=[1]^\perp \cap W_1$ and then one simply has ${\Vert P_{Z_1} \mu\Vert}^2=\sum(\mu_i - \bar \mu)^2$.

Now do you define the population $R^2$ of the submodel $H_0$ as the asymptotic value of the $R^2$ calculated with respect to model $H_0$ but under the distributional assumption of model $H_1$ ? The asymptotic value (if there is one) seems more difficult to find.

$\endgroup$
  • $\begingroup$ Thanks Stéphane. I'll have to have a think about what you are saying. With regards to your question. I assume that the true data generating process is not known but that it is the same fo both models, but that there is a true proportion of variance explained by the linear regression in model 1 and model 2. $\endgroup$ – Jeromy Anglim Jun 30 '13 at 7:48
  • $\begingroup$ @JeromyAnglim Formula (A3) of this paper is a particular case of my formula for the one-way ANOVA model. So my formula should be the general definition of the population $R^2$, but this is not what you are using in your OP. $\endgroup$ – Stéphane Laurent Jun 30 '13 at 21:37
  • 1
    $\begingroup$ @JeromyAnglim The study of this paper seems to be close to what you're looking for (with random predictors). $\endgroup$ – Stéphane Laurent Jun 30 '13 at 21:48
  • $\begingroup$ Thanks. The Algina, Keselman and Penfield paper looks very useful. I added some comments to my answer about it. $\endgroup$ – Jeromy Anglim Jul 1 '13 at 7:21
  • $\begingroup$ @JeromyAnglim So what is the assumption about the predictors ? They are generated according to a multivariate Gaussian distribution ? $\endgroup$ – Stéphane Laurent Jul 3 '13 at 6:12
1
$\begingroup$

Rather than answer the question you asked, I'm going to ask why you ask that question. I assume you want to know whether

mod.small <- lm(y ~ x1a + x1b + x1c, data=x)

is at least as good as

mod.large <- lm(y ~ ., data=x)

at explaining y. Since these models are nested, the obvious way to answer this question would seem to be to run an analysis of variance comparing them, in the same way as you might run an analysis of deviance for two GLMs, like

anova(mod.small, mod.large)

Then you could use the sample R-square improvement between models as your best guess at what the fit improvement would be in the population, always assuming you can make sense of population R-squared. Personally I'm not sure I can, but with this it doesn't matter either way.

More generally, if you're interested in population quantities you're presumably interested in generalisation so a sample fit measure is not quite what you want, however 'corrected'. For example, cross-validation of some quantity that estimates the sort and quantity of actual errors you could expect to make out of sample, like MSE, would seem to get at what you want.

But it's quite possible I'm missing something here...

$\endgroup$
  • $\begingroup$ I appreciate your answer, and it may well be good advice for others. But my research context means that I am legitimately interested in delta-rho square. While most statisticians are often more concerned with the predictive utility of a model (e.g., cross-validated delta r-square), I am a psychological scientist and am specifically interested in the population property. Furthermore, I am not interested in statistical significance of the improvement. I am interested in the size of the improvement. And I find that delta-r-square is a useful metric for indexing that size of improvement. $\endgroup$ – Jeromy Anglim Jul 1 '13 at 7:15
  • $\begingroup$ With regards to MSE, different studies in psychology use measures on very different metrics. Thus, there is an attraction, right or wrong, to standardised measures such as r-square. $\endgroup$ – Jeromy Anglim Jul 1 '13 at 7:24
  • $\begingroup$ Fair enough, particularly on MSE. I remain a bit confused by the interest in bootstrapping and population inference but the lack of interest in testing since, perhaps naively, these seem to be equivalent concerns differently addressed. I'm also having difficulty tightly distinguishing out of sample prediction from inference to a population, but that's probably pre-coffee knee-jerk bayesianism (where prediction is just another population inference problem) getting in the way. $\endgroup$ – conjugateprior Jul 1 '13 at 8:59
  • $\begingroup$ Perhaps I spoke a bit quickly. In my research context, there is often plenty of evidence that the delta-rho-square is greater than zero. The question of interest is what is the degree of increase. I.e., it is a trivial increase or a theoretically meaningful increase. Thus, confidence or credible intervals give me an estimate of the uncertainty around that increase. I haven't yet reconciled what I'm doing here with my understanding of Bayesian statistics, but I'd like to. $\endgroup$ – Jeromy Anglim Jul 1 '13 at 9:19
1
$\begingroup$

The following represent a few possibilities for calculating confidence intervals on $\rho^2$.

Double adjusted r-square bootstrap

My current best guess at an answer is to do a double adjusted r-square bootstrap. I've implemented the technique. It involves the following:

  • Generate a set of bootstrap samples from the current data.
  • For each bootstrapped sample:
    • calculate first adjusted r-square for the two models
    • calculate second adjusted r-square on the adjusted r-square values from the previous step
    • Subtract model2 from model1 second adjusted r-square values to get an estimate of $\Delta \rho^2$.

The rationale is that the first adjusted r-square removes the bias introduced by bootsrapping (i.e., bootstrapping assumes that the sample r-square is the population r-square). The second adjusted r-square performs the standard correction that is applied to a normal sample to estimate population r-square.

At this point, all I can see is that applying this algorithm generates estimates that seem about right (i.e., the mean theta_hat in the bootstrap is very close to the sample theta_hat). The standard error aligns with my intuition. I haven't yet tested whether it provides proper frequentist coverage where the data generating process is known, and I'm also not entirely sure at this point how the argument could be justified from first principles

If anyone sees any reasons why this approach would be problematic, I'd be grateful to hear about it.

Simulation by Algina et al

Stéphane mentioned the article by Algina, Keselman and Penfield. They performed a simulation study to examine the 95% confidence interval coverage of bootstrapping and asymptotic methods for estimating $\Delta \rho^2$. Their bootstrapping methods involved only a single application of adjusted r-square, rather than the double adjustment of r-square that I mention above. They found that bootstrap estimates only provided good coverage when the number of additional predictors in the full model was one or perhaps two. It is my hypothesis that this is because as the number of predictors increases, so would the difference between the single and double adjusted r-square bootstrap.

Smithson (2001) on using the noncentrality parameter

Smithson (2001) discusses calculating confidence intervals for the partial $R^2$ based on the non-centrality parameter. See pages 615 and 616 in particular. He suggests that "it is straightforward to construct a CI for $f^2$ and partial $R^2$ but not for the squared semipartial correlation." (p.615)

References

  • Algina, J., Keselman, H. J., & Penfield, R. D. Confidence Intervals for the Squared Multiple Semipartial Correlation Coefficient. PDF
  • Smithson, M. (2001). Correct confidence intervals for various regression effect sizes and parameters: The importance of noncentral distributions in computing intervals. Educational and Psychological Measurement, 61(4), 605-632.
$\endgroup$
  • 1
    $\begingroup$ It seems that nobody here (including you) knows the definition of your population R-squared. Hence IMHO this is a seriously problematic approach. $\endgroup$ – Stéphane Laurent Jul 1 '13 at 5:35
  • $\begingroup$ @StéphaneLaurent Thanks for that. I confess that up to this point I have not seen population r-square as a property of contention. For example, I could propose a data generating process and there would be an r-square that is approached as my simulation sample size approaches infinity. And likewise I assume that there is a data generating process for my data, and therefore if it were possible to get an infinite sample, I could calculate the true population r-square. $\endgroup$ – Jeromy Anglim Jul 1 '13 at 7:20
  • $\begingroup$ Yes but I'm under the impression you also assume a generating process for the predictors. I cannot figure how this could make sense for a general linear model. $\endgroup$ – Stéphane Laurent Jul 1 '13 at 8:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.