9
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Could you clarify why you can't include both predictors in the model? In your case, X1 and X2 would probably be correlated, as the brown spots probably increase with increasing time lying on the table. $\endgroup$ – COOLSerdash Jun 25 '13 at 8:45
  • $\begingroup$ I'm interesting in testing whether X1 or X2 is the better measurement method. If including them both in one model can answer that question, there is no problem doing that. Obviously they are both correlated as they measure the same thing. $\endgroup$ – Rodin Jun 25 '13 at 8:57
  • $\begingroup$ I would like to say: when trying to measure banana tastiness, measuring how long it's been lying on the table is a better way to determine this than counting the number of brown spots (p < 0.05). $\endgroup$ – Rodin Jun 25 '13 at 9:04
7
$\begingroup$

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)

SPLOM

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ This is a practical detailed answer explaining well what could be done. What would have been even more striking would have been an example in which scatter plots, $R^2$ and $AIC_c$ lead to different answers. $\endgroup$ – Nick Cox Jun 25 '13 at 13:02
  • $\begingroup$ @NickCox That would indeed be interesting! If the models are linear and of equal complexity, and the errors have a Gaussian distribution, $R^2, AIC$ and $BIC$ should give the same answers. The truth is that I don't know how to produce an example ad-hoc where scatterplots, $R^2$ and $AIC_{c}$ lead to different answers (maybe someone could help me?). Sorry. $\endgroup$ – COOLSerdash Jun 25 '13 at 13:27
  • $\begingroup$ You need not apologise; I see it as good news whenever different sensible methods imply the same conclusion! $\endgroup$ – Nick Cox Jun 25 '13 at 13:29
  • $\begingroup$ Excellent. The Cox test is exactly what I wanted. Unfortunately, my models are linear-mixed-effect models fitted with the lme4 package, which are not directly supported by the lmtest package. Before I dive into the literature written about cox-like tests with LMEs, does anyone know a readily available R-package to do it? $\endgroup$ – Rodin Jun 25 '13 at 14:49
  • $\begingroup$ @Rodin No, I don't know any R package that could do that. Maybe this post can give you further guidance. $\endgroup$ – COOLSerdash Jun 25 '13 at 14:52
3
$\begingroup$

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?".

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ +1, good points. I have explicitly stated that one should not use significance testing when dealing with the $AIC$ and alike. $\endgroup$ – COOLSerdash Jun 25 '13 at 11:28
  • 1
    $\begingroup$ Always good to take a step back and remember this so +1 but I am wondering if this might in fact just be the odd case where this is not particularly useful. Is it really quite likely that one model will be clearly better than the other when the two tentative predictors are different measures of the same thing as opposed to substantively different variables? Leaving aside the bananas for a minute, think slightly different questionnaires or a ruler vs. a laser rangefinder. One measure has to be extremely deficient for non-linearity to show up in one case but not the other. $\endgroup$ – Gala Jun 25 '13 at 11:53
  • 2
    $\begingroup$ I agree, but I feel happy to extend my answer to residual plots, etc. More generally, any way to judge this could run into the same objection, e.g. significance tests might fail to reject equal performance, $R^2$ might be practically identical. Having published on comparison of different measurement methods with real data, I do assert that choice of method will usually find a best method, but if two methods appear about equally good, then that's the answer, not an insoluble problem leading to methodological anguish. $\endgroup$ – Nick Cox Jun 25 '13 at 12:02
  • $\begingroup$ I agree that the best way is to clearly present the data and the fits by the models. This way the reader can decide for themselves whether to accept any statement about one being better or not. However, I'm afraid reviewers will demand a significance test, purely out of a knee-jerk reaction. Gaël's comment about night and day is not so far off. $\endgroup$ – Rodin Jun 25 '13 at 14:26
  • 1
    $\begingroup$ That was me on night and day.... $\endgroup$ – Nick Cox Jun 25 '13 at 14:33
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$

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.