Linear model: comparing predictive power of two different measurement methods

Question

I'm interested in predicting Y and am studying different two measurement techniques X1 and X2. It could be for instance that I want to predict the tastiness of a banana, either by measuring how long it has been lying on the table, or by measuring the number of brown spots on the banana.

I want to know which one of the measuring techniques is better, should I choose to perform only one.

I can create a linear model in R:

m1 = lm(Y ~ X1)
m2 = lm(Y ~ X2)

Now let's say X1 is a superior predictor of banana tastiness than X2. When calculating the $R^2$ of the two models, the $R^2$ of model m1 is clearly higher than model m2. Before writing a paper on how method X1 is better than X2, I want to have some sort of indication that the difference is not by chance, possibly in the form of a p-value.

How would one go about this? How to do it when I'm using different brands of bananas and move to a Linear Mixed Effect model that incoporates banana brand as a random effect?

Answer 1

Later

One thing I want to add after hearing that you have linear mixed effect models: The $AIC, AIC_{c}$ and $BIC$ can still be used to compare the models. See this paper, for example. From other similar questions on the site, it seems that this paper is crucial.

---

Original answer

What you basically want is to compare two non-nested models. Burnham and Anderson Model selection and multimodel inference discuss this and recommend using the $AIC$, $AIC_{c}$ or $BIC$ etc. as the traditional likelihood ratio test is only applicable in nested models. They explicitly state that the information-theoretic criteria such as the $AIC, AIC_{c}, BIC$ etc. are not tests and that the word "significant" should be avoided when reporting the results.

Based on this and this answers, I recommend these approaches:

1. Make a scatterplot matrix (SPLOM) of your dataset including smoothers: pairs(Y~X1+X2, panel = panel.smooth, lwd = 2, cex = 1.5, col = "steelblue", pch=16). Check if the lines (the smoothers) are compatible with a linear relationship. Refine the model if necessary.
2. Compute the models m1 and m2. Do some model checks (residuals etc.): plot(m1) and plot(m2).
3. Compute the $AIC_{c}$ ($AIC$ corrected for small sample sizes) for both models and calculate the absolute difference between the two $AIC_{c}$s. The R package pscl provides the function AICc for this: abs(AICc(m1)-AICc(m2)). If this absolute difference is smaller than 2, the two models are basically indistinguishable. Otherwise prefer the model with the lower $AIC_{c}$.
4. Compute likelihood ratio tests for non-nested models. The R package lmtest has the functions coxtest (Cox test), jtest (Davidson-MacKinnon J test) and encomptest (encompassing test of Davidson & MacKinnon).

Some thoughts: If the two banana-measures are really measure the same thing, they both may be equally suited for prediction and there might not be a "best" model.

This paper might also be helpful.

Here is an exmple in R:

#==============================================================================
# Generate correlated variables
#==============================================================================

set.seed(123)

R <- matrix(cbind(
  1   , 0.8 , 0.2,
  0.8 , 1   , 0.4,
  0.2 , 0.4 , 1),nrow=3) # correlation matrix
U <- t(chol(R))
nvars <- dim(U)[1]
numobs <- 500
set.seed(1)
random.normal <- matrix(rnorm(nvars*numobs,0,1), nrow=nvars, ncol=numobs);
X <- U %*% random.normal
newX <- t(X)
raw <- as.data.frame(newX)
names(raw) <- c("response","predictor1","predictor2")

#==============================================================================
# Check the graphic
#==============================================================================

par(bg="white", cex=1.2)
pairs(response~predictor1+predictor2, data=raw, panel = panel.smooth,
      lwd = 2, cex = 1.5, col = "steelblue", pch=16, las=1)

The smoothers confirm the linear relationships. This was intended, of course.

#==============================================================================
# Calculate the regression models and AICcs
#==============================================================================

library(pscl)

m1 <- lm(response~predictor1, data=raw)
m2 <- lm(response~predictor2, data=raw)

summary(m1)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.004332   0.027292  -0.159    0.874    
predictor1   0.820150   0.026677  30.743   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6102 on 498 degrees of freedom
Multiple R-squared:  0.6549,    Adjusted R-squared:  0.6542 
F-statistic: 945.2 on 1 and 498 DF,  p-value: < 2.2e-16

summary(m2)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.01650    0.04567  -0.361    0.718    
predictor2   0.18282    0.04406   4.150 3.91e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.021 on 498 degrees of freedom
Multiple R-squared:  0.03342,   Adjusted R-squared:  0.03148 
F-statistic: 17.22 on 1 and 498 DF,  p-value: 3.913e-05

AICc(m1)
[1] 928.9961

AICc(m2)
[1] 1443.994

abs(AICc(m1)-AICc(m2))
[1] 514.9977

#==============================================================================
# Calculate the Cox test and Davidson-MacKinnon J test
#==============================================================================

library(lmtest)

coxtest(m1, m2)

Cox test

Model 1: response ~ predictor1
Model 2: response ~ predictor2
                Estimate Std. Error   z value  Pr(>|z|)    
fitted(M1) ~ M2   17.102     4.1890    4.0826 4.454e-05 ***
fitted(M2) ~ M1 -264.753     1.4368 -184.2652 < 2.2e-16 ***

jtest(m1, m2)

J test

Model 1: response ~ predictor1
Model 2: response ~ predictor2
                Estimate Std. Error t value  Pr(>|t|)    
M1 + fitted(M2)  -0.8298   0.151702  -5.470 7.143e-08 ***
M2 + fitted(M1)   1.0723   0.034271  31.288 < 2.2e-16 ***

The $AIC_{c}$ of the first model m1 is clearly lower and the $R^{2}$ is much higher.

Important: In linear models of equal complexity and Gaussian error distribution, $R^2, AIC$ and $BIC$ should give the same answers (see this post). In nonlinear models, the use of $R^2$ for model performance (goodness of fit) and model selection should be avoided: see this post and this paper, for example.

Answer 2

There's a good 19th century answer in risk of neglect here. To compare two different straight-line fits, plot the data and the fitted lines and think about what you see. It's quite likely that one model will be clearly better, and that need not mean higher $R^2$. For example, it is possible that a straight-line model is qualitatively wrong in one or other case. Even better, the data and fit may suggest a better model. If the two models appear about equally good or poor, that's another answer.

The banana example is presumably facetious here, but I would not expect straight-line fits to work well at all....

The inferential machinery wheeled out by others in answer is a thing of intellectual beauty, but sometimes you don't need a state-of-the-art sledgehammer to crack a nut. Sometimes it seems that anyone publishing that night is darker than day would always have some one asking "Have you tested that formally? What is your P-value?".

Answer 3

Do a Cox test for non-nested models.

y <- rnorm( 10 )
x1 <- y + rnorm( 10 ) / 2
x2 <- y + rnorm( 10 )

lm1 <- lm( y ~ x1 )
lm2 <- lm( y ~ x2 )

library( lmtest )

coxtest( lm1, lm2 )
?coxtest

(you will find references to other test).

See also this comment and this question. In especially, consider using AIC/BIC.