Quelle est la puissance du test de régression F?

11

Le test F classique pour des sous-ensembles de variables en régression multilinéaire a la forme où est la somme des erreurs quadratiques sous le modèle "réduit", qui nichent à l'intérieur du "grand" modèle , et sont les degrés de liberté du deux modèles. Dans l'hypothèse nulle selon laquelle les variables supplémentaires du «grand» modèle n'ont pas de pouvoir explicatif linéaire, la statistique est distribuée sous la forme d'un F avec et df_B degrés de liberté.

F=(SSE(R)SSE(B))/(dfRdfB)SSE(B)/dfB,
B d f d f R - d f B d f BSSE(R)BdfdfRdfBdfB

Quelle est cependant la répartition sous l'alternative? Je suppose que c'est un F non central (j'espère pas doublement non central), mais je ne trouve aucune référence sur ce qu'est exactement le paramètre de non-centralité. Je vais deviner que cela dépend des vrais coefficients de régression β , et probablement de la matrice de conception X , mais au-delà, je ne suis pas si sûr.

shabbychef
la source

Réponses:

9

Le paramètre de non-centralité est δ2 , la projection pour le modèle restreint est Pr , β est le vecteur de vrais paramètres, X est la matrice de conception pour le modèle sans restriction (vrai), ||x||est la norme:

δ2=||XβPrXβ||2σ2

Vous pouvez lire la formule comme ceci: est le vecteur des valeurs attendues sous condition de la matrice de conception . Si vous traitez comme un vecteur de données empirique , alors sa projection sur le sous-espace du modèle restreint est , ce qui vous donne la prédiction du modèle restreint pour ces "données". Par conséquent, est analogue à et vous donne l'erreur de cette prédiction. Par conséquent donne la somme des carrés de cette erreur. Si le modèle restreint est vrai, alorsX X β y P r X β y X β - P r X β y - y | | X β - P r X β | | 2 X β X rE(y|X)=XβXXβyPrXβy^XβPrXβyy^||XβPrXβ||2Xβse trouve déjà dans le sous-espace défini par et , de sorte que le paramètre de non-centralité est .Xr0PrXβ=Xβ0

Vous devriez trouver cela à Mardia, Kent & Bibby. (1980). Analyse multivariée.

caracal
la source
génial! faut-il mettre la norme au carré? Sinon, il semble que les unités comptent? Vous dites que c'est "la somme des carrés", donc je pense que c'est la norme au carré ..
shabbychef
@shabbychef Bien sûr, vous avez raison, merci d'avoir attrapé ça!
caracal
7

J'ai confirmé la réponse de @ caracal avec une expérience de Monte Carlo. J'ai généré des instances aléatoires à partir d'un modèle linéaire (avec la taille aléatoire), calculé la statistique F et calculé la valeur p en utilisant le paramètre de non-centralité puis j'ai tracé le cdf empirique de ces valeurs de p. Si le paramètre de non-centralité (et le code!) Est correct, je devrais obtenir un cdf presque uniforme, ce qui est le cas:

δ2=||Xβ1Xβ2||2σ2,

CDF empirique de ce qui devrait être normal

Voici le code R (pardonnez le style, j'apprends encore):

#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))
shabbychef
la source
2
+1 pour le suivi avec le code. C'est toujours bon de voir ça.
mpiktas