¿Cuál es el poder de la regresión de la prueba de F?

Question

La clásica prueba F de subconjuntos de variables en la regresión multilineal tiene la forma $$ F = \frac{(\mbox{ESS}(R) - \mbox{ESS}(B))/(df_R - df_B)}{\mbox{ESS}(B)/df_B}, $$ donde $\mbox{SSE}(R)$ es la suma de los cuadrados de los errores bajo el "reducido" del modelo, que anidan en el interior de la "gran" modelo de $B$, e $df$ son los grados de libertad de los dos modelos. Bajo la hipótesis nula de que las variables extra en la 'gran' modelo no lineal poder explicativo, el estadístico se distribuye como una F con $df_R - df_B$ $df_B$ grados de libertad.

¿Cuál es la distribución, sin embargo, bajo la alternativa? Supongo que es un no-central F (espero que no sea doblemente no central), pero no puedo encontrar ninguna referencia sobre lo que es exactamente el parámetro no-centralidad. Yo voy a adivinar que depende de la verdadera coeficientes de regresión $\beta$, y probablemente en el diseño de la matriz de $X$, pero más allá de eso no estoy tan seguro.

Answer 1

El noncentrality parámetro es $\delta^{2}$, la proyección para el modelo restringido es $P_{r}$, $\beta$ es el vector de los verdaderos parámetros, $X$ es el diseño de la matriz de restricciones (true) modelo, $|| x ||$ es la norma:

$$ \delta^{2} = \frac{|| X \beta - P_{r} X \beta ||^{2}}{\sigma^{2}} $$

Usted puede leer la fórmula como esta: $E(y | X) = X \beta$ es el vector de valores esperados condicionales en el diseño de la matriz de $X$. Si usted trata a $X \beta$ como empírico en el vector de datos $y$, luego de su proyección en el modelo restringido subespacio es $P_{r} X \beta$, lo que le da la predicción de $\hat{y}$ de la restringida modelo para que los "datos". En consecuencia, $X \beta - P_{r} X \beta$ es análoga a $y - \hat{y}$ y te da el error de la predicción. Por lo tanto $|| X \beta - P_{r} X \beta ||^{2}$ da la suma de los cuadrados del error. Si el modelo restringido es cierto, entonces, $X \beta$ ya está en el subespacio definido por $X_{r}$, e $P_{r} X \beta = X \beta$, de tal manera que el noncentrality parámetro es $0$.

Usted debe encontrar esto en Mardia, Kent & Bibby. (1980). En El Análisis Multivariante.

Answer 2

He confirmado @caracal la respuesta con un experimento de Monte Carlo. Yo aleatorio generado instancias de un modelo lineal (con el tamaño aleatorio), se calcula el estadístico F y calcula el p-valor mediante el parámetro no-centralidad de $$ \delta^2 = \frac{||X\beta_1 - X\beta_2||^2}{\sigma^2}, $$ A continuación, he trazado la cdf empírica de estos p-valores. Si el parámetro no-centralidad (y el código!) es correcto, yo debería obtener una cerca de uniforme cdf, que es el caso de:

Aquí está el código R (perdón por el estilo, todavía estoy aprendiendo):

#sum of squares
sum2 <- function(x) { return(sum(x * x)) }
#random integer between n and 2n
rint <- function(n) { return(ceiling(runif(1,min=n,max=2*n))) }
#generate random instance from linear model plus noise.
#n observations of p2 vector
#regress against all variables and against a subset of p1 of them
#compute the F-statistic for the test of the p2-p1 marginal variables
#compute the p-value under the putative non-centrality parameter
gend <- function(n,p1,p2,sig = 1) {
 beta2 <- matrix(rnorm(p2,sd=0.1),nrow=p2)
 beta1 <- matrix(beta2[1:p1],nrow=p1)
 X <- matrix(rnorm(n*p2),nrow=n,ncol=p2)
 yt1 <- X[,1:p1] %*% beta1
 yt2 <- X %*% beta2
 y <- yt2 + matrix(rnorm(n,mean=0,sd=sig),nrow=n)
 ncp <- (sum2(yt2 - yt1)) / (sig ** 2)
 bhat2 <- lm(y ~ X - 1)
 bhat1 <- lm(y ~ X[,1:p1] - 1)
 SSE1 <- sum2(bhat1$residual)
 SSE2 <- sum2(bhat2$residual)
 df1 <- bhat1$df.residual
 df2 <- bhat2$df.residual
 Fstat <- ((SSE1 - SSE2) / (df1 - df2)) / (SSE2 / bhat2$df.residual)
 pval <- pf(Fstat,df=df1-df2,df2=df2,ncp=ncp)
 return(pval)
}
#call the above function, but randomize the problem size (within reason)
genr <- function(n,p1,p2,sig=1) {
 use.p1 <- rint(p1)
 use.p2 <- use.p1 + rint(p2 - p1)
 return(gend(n=rint(n),p1=use.p1,p2=use.p2,sig=sig+runif(1)))
}
ntrial <- 4096
ssize <- 256
z <- replicate(ntrial,genr(ssize,p1=4,p2=10))
plot(ecdf(z))