Modelo lineal: comparación del poder predictivo de dos métodos de medición diferentes

Estoy interesado en predecir Yy estoy estudiando diferentes técnicas de medición X1y X2. Podría ser, por ejemplo, que quiero predecir el sabor de un plátano, ya sea midiendo cuánto tiempo ha estado sobre la mesa o midiendo la cantidad de manchas marrones en el plátano.

Quiero saber cuál de las técnicas de medición es mejor, si elijo realizar solo una.

Puedo crear un modelo lineal en R:

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

Ahora digamos que X1es un predictor superior de sabor a plátano que X2. Al calcular el de los dos modelos, el del modelo es claramente más alto que el modelo . Antes de escribir un artículo sobre cómo es mejor el método , quiero tener algún tipo de indicación de que la diferencia no es casual, posiblemente en forma de un valor p.$R^2$$R^2$m1m2X1X2

¿Cómo se podría hacer esto? ¿Cómo hacerlo cuando estoy usando diferentes marcas de plátanos y pasar a un modelo de efecto mixto lineal que incorpora la marca de plátanos como un efecto aleatorio?

Más tarde

Una cosa que quiero agregar después de escuchar que tiene modelos lineales de efectos mixtos: el y todavía se pueden usar para comparar los modelos. Ver este artículo , por ejemplo. De otras preguntas similares en el sitio, parece que este documento es crucial.$AIC, AIC_{c}$$BIC$

Respuesta original

Lo que básicamente quiere es comparar dos modelos no anidados. La selección de modelos de Burnham y Anderson y la inferencia multimodelo discuten esto y recomiendan el uso de , o etc., ya que la prueba de razón de probabilidad tradicional solo es aplicable en modelos anidados. Indican explícitamente que los criterios teóricos de la información, como el etc. no son pruebas y que la palabra "significativo" debe evitarse al informar los resultados.$AIC$$AIC_{c}$$BIC$$AIC, AIC_{c}, BIC$

Basado en esto y en estas respuestas, recomiendo estos enfoques:

1. Hacer un diagrama de dispersión matricial (SPLOM) del conjunto de datos que incluye alisadores: pairs(Y~X1+X2, panel = panel.smooth, lwd = 2, cex = 1.5, col = "steelblue", pch=16). Compruebe si las líneas (los suavizadores) son compatibles con una relación lineal. Refina el modelo si es necesario.
2. Calcular los modelos m1y m2. Haga algunas comprobaciones del modelo (residuos, etc.): plot(m1)y plot(m2).
3. Calcule el ( corregido para tamaños de muestra pequeños) para ambos modelos y calcule la diferencia absoluta entre los dos s. El paquete proporciona la función para esto: . Si esta diferencia absoluta es menor que 2, los dos modelos son básicamente indistinguibles. De lo contrario prefiere el modelo con el menor .$AIC_{c}$$AIC$$AIC_{c}$R psclAICcabs(AICc(m1)-AICc(m2))$AIC_{c}$
4. Calcule las pruebas de razón de probabilidad para modelos no anidados. El R paquetelmtest tiene las funciones coxtest(prueba de Cox), jtest( prueba de Davidson-MacKinnon J) y encomptest(prueba de Davidson y MacKinnon).

Algunas reflexiones: si las dos medidas de banano realmente miden lo mismo, ambas pueden ser igualmente adecuadas para la predicción y puede que no haya un "mejor" modelo.

Este documento también puede ser útil.

Aquí hay un ejemplo en 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)

Los suavizadores confirman las relaciones lineales. Esto estaba destinado, por supuesto.

#==============================================================================
# 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 ***

El del primer modelo es claramente más bajo y el es mucho más alto.$AIC_{c}$m1$R^{2}$

Importante: en modelos lineales de igual complejidad y distribución de errores gaussianos , y deberían dar las mismas respuestas (ver esta publicación ). En los modelos no lineales , se debe evitar el uso de para el rendimiento del modelo (bondad de ajuste) y la selección del modelo : consulte esta publicación y este documento , por ejemplo.$R^2, AIC$$BIC$$R^2$

Aquí hay una buena respuesta del siglo XIX en riesgo de negligencia. Para comparar dos ajustes de línea recta diferentes, trace los datos y las líneas ajustadas y piense en lo que ve. Es muy probable que un modelo sea claramente mejor, y eso no necesariamente significa más alto$R^2$. Por ejemplo, es posible que un modelo lineal sea cualitativamente incorrecto en uno u otro caso. Aún mejor, los datos y el ajuste pueden sugerir un mejor modelo. Si los dos modelos parecen igualmente buenos o malos, esa es otra respuesta.

El ejemplo del plátano es presumiblemente gracioso aquí, pero no esperaría que los ajustes en línea recta funcionen bien ...

La maquinaria inferencial extraída por otros en respuesta es una cosa de belleza intelectual, pero a veces no se necesita un mazo de última generación para romper una nuez. A veces parece que cualquiera que publique esa noche es más oscuro que el día siempre tendría a alguien preguntando "¿Lo ha probado formalmente? ¿Cuál es su valor P?".

Haga una prueba de Cox para modelos no anidados.

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

(encontrará referencias a otras pruebas).

Vea también este comentario y esta pregunta . En especial, considere usar AIC / BIC.