Tengo algunos datos que están altamente correlacionados. Si ejecuto una regresión lineal, obtengo una recta de regresión con una pendiente cercana a uno (= 0,93). Lo que me gustaría hacer es comprobar si esta pendiente es significativamente diferente de 1,0. Mi expectativa es que no lo es. En otras palabras, me gustaría cambiar la hipótesis nula de la regresión lineal de una pendiente cero a una pendiente uno. ¿Es éste un enfoque sensato? También te agradecería mucho que incluyeras algún código de R en tu respuesta para poder implementar este método (¡o uno mejor que sugieras!). Gracias.
Respuestas
¿Demasiados anuncios?set.seed(20); y = rnorm(20); x = y + rnorm(20, 0, 0.2) # generate correlated data
summary(lm(y ~ x)) # original model
summary(lm(y ~ x, offset= 1.00*x)) # testing against slope=1
summary(lm(y-x ~ x)) # testing against slope=1
Salidas:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01532 0.04728 0.324 0.75
x 0.91424 0.04128 22.148 1.64e-14 ***
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01532 0.04728 0.324 0.7497
x -0.08576 0.04128 -2.078 0.0523 .
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.01532 0.04728 0.324 0.7497
x -0.08576 0.04128 -2.078 0.0523 .
Su hipótesis puede expresarse como $R\beta=r$ donde $\beta$ son sus coeficientes de regresión y $R$ es la matriz de restricción con $r$ las restricciones. Si nuestro modelo es
$$y=\beta_0+\beta_1x+u$$
entonces para la hipótesis $\beta_1=0$ , $R=[0,1]$ y $r=1$ .
Para este tipo de hipótesis puede utilizar linearHypothesis
del paquete coche :
set.seed(20); y = rnorm(20); x = y + rnorm(20, 0, 0.2) # generate correlated data
mod <- lm(y ~ x)) # original model
> linearHypothesis(mod,matrix(c(0,1),nrow=1),rhs=c(1))
Linear hypothesis test
Hypothesis:
x = 1
Model 1: restricted model
Model 2: y ~ x
Res.Df RSS Df Sum of Sq F Pr(>F)
1 19 0.96022
2 18 0.77450 1 0.18572 4.3162 0.05234 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Parece que sigues intentando rechazar una hipótesis nula. Hay un montón de problemas con eso, no menos importante de los cuales es que es posible que usted no tiene suficiente poder para ver que es diferente de 1. Suena como que no le importa que la pendiente es 0,07 diferente de 1. Pero ¿y si realmente no se puede saber? ¿Qué pasa si en realidad está estimando una pendiente que varía salvajemente y en realidad puede estar muy lejos de 1 con algo así como un intervalo de confianza de ±0,4? Su mejor táctica aquí no es cambiar la hipótesis nula, sino realmente hablar razonablemente sobre una estimación de intervalo. Si aplica el comando confint() a su modelo puede obtener un intervalo de confianza del 95% alrededor de su pendiente. Entonces puede utilizarlo para discutir la pendiente que obtuvo. Si 1 está dentro del intervalo de confianza, puede afirmar que está dentro del rango de valores que cree probable que contenga el valor verdadero. Pero lo más importante es que también puedes decir cuál es ese intervalo de valores.
El objetivo de las pruebas es rechazar la hipótesis nula, no confirmarla. El hecho de que no haya una diferencia significativa, no es en absoluto una prueba de la ausencia de una diferencia significativa. Para ello, tendrá que definir qué tamaño del efecto considera razonable para rechazar la nula.
Comprobar si la pendiente es significativamente diferente de 1 no es tan difícil, basta con comprobar si la diferencia $slope - 1$ difiere significativamente de cero. A mano sería algo como :
set.seed(20); y = rnorm(20); x = y + rnorm(20, 0, 0.2)
model <- lm(y~x)
coefx <- coef(summary(model))[2,1]
seslope <- coef(summary(model))[2,2]
DF <- model$df.residual
# normal test
p <- (1 - pt(coefx/seslope,DF) )*2
# test whether different from 1
p2 <- (1 - pt(abs(coefx-1)/seslope,DF) )*2
Ahora debe ser consciente del hecho de que el tamaño del efecto para el que una diferencia se convierte en significativa, es
> qt(0.975,DF)*seslope
[1] 0.08672358
siempre que tengamos un estimador decente del error estándar de la pendiente. Por lo tanto, si decide que una diferencia significativa sólo debe detectarse a partir de 0,1, puede calcular la DF necesaria de la siguiente manera :
optimize(
function(x)abs(qt(0.975,x)*seslope - 0.1),
interval=c(5,500)
)
$minimum
[1] 6.2593
Eso sí, esto depende bastante de la estimación del seslope. Para obtener una mejor estimación de la seslope, podrías hacer un remuestreo de tus datos. Una forma ingenua sería :
n <- length(y)
seslope2 <-
mean(
replicate(n,{
id <- sample(seq.int(n),1)
model <- lm(y[-id]~x[-id])
coef(summary(model))[2,2]
})
)
poniendo seslope2 en la función de optimización, devuelve :
$minimum
[1] 6.954609
Todo esto le dirá que su conjunto de datos devolverá un resultado significativo más rápido de lo que usted considera necesario, y que sólo necesita 7 grados de libertad (en este caso 9 observaciones) si quiere estar seguro de que no significativo significa lo que usted quiere que signifique.
Simplemente no se pueden hacer afirmaciones de probabilidad o verosimilitud sobre el parámetro utilizando un intervalo de confianza, esto es un paradigma bayesiano.
Lo que John dice es confuso porque hay una equivalencia entre CIs y Pvalues, así que a un 5%, decir que tu CI incluye 1 es equivalente a decir que Pval>0.05.
linearHypothesis permite probar restricciones diferentes de la beta estándar=0