Pourquoi lm et biglm dans R donnent-ils des valeurs de p différentes pour les mêmes données?

Voici un petit exemple:

MyDf<-data.frame(x=c(1,2,3,4), y=c(1.2, .7, -.5, -3))

Maintenant avec le base::lm:

> lm(y~x, data=MyDf) %>% summary

Call:
lm(formula = y ~ x, data = MyDf)

Residuals:
    1     2     3     4 
-0.47  0.41  0.59 -0.53 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   3.0500     0.8738   3.491   0.0732 .
x            -1.3800     0.3191  -4.325   0.0495 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7134 on 2 degrees of freedom
Multiple R-squared:  0.9034,    Adjusted R-squared:  0.8551 
F-statistic: 18.71 on 1 and 2 DF,  p-value: 0.04952

Maintenant, essayez la même chose avec biglmle biglmpackage:

XX<-biglm(y~x, data=MyDf) 
print(summary(XX), digits=5)

Large data regression model: biglm(y ~ x, data = MyDf)
Sample size =  4 
             Coef     (95%      CI)      SE       p
(Intercept)  3.05  1.30243  4.79757 0.87378 0.00048
x           -1.38 -2.01812 -0.74188 0.31906 0.00002

Notez que nous avons besoin de printet digitspour voir la valeur de p. Les coefficients et les erreurs standard sont les mêmes, mais les valeurs de p sont très différentes. Pourquoi cela est-il ainsi?

r regression p-value linear-model John Paul
la source

+1 indice: comparer pt(-3.491, 2)*2à pnorm(-3.491)*2, par exemple.

whuber

@whuber Merci. Il s'agit donc essentiellement d'un problème de distribution t par rapport à la distribution normale. L'idée est-elle que la distribution normale a plus de sens pour les grands ensembles de données typiques de biglm?

John Paul

Je pense que l'idée est que la normale n'est pas si différente de t avec une valeur élevée. Essayez l'exemple du premier commentaire, mais changez pt (-3.491, 2) * 2 en pt (-3.491, 2e3) * 2.

ν

$\nu$

Andrey Kolyadin

Pour voir quelles valeurs de p sont correctes (le cas échéant), répétons le calcul pour les données simulées dans lesquelles l'hypothèse nulle est vraie. Dans le cadre actuel, le calcul est un ajustement des moindres carrés aux données (x, y) et l'hypothèse nulle est que la pente est nulle. Dans la question, il y a quatre valeurs x 1,2,3,4 et l'erreur estimée est d'environ 0,7, alors incorporons cela dans la simulation.

Voici la configuration, écrite pour être compréhensible par tout le monde, même ceux qui ne le connaissent pas R.

beta <- c(intercept=0, slope=0)
sigma <- 0.7
x <- 1:4
y.expected <-  beta["intercept"] + beta["slope"] * x

La simulation génère des erreurs indépendantes, les ajoute à y.expected, appelle lmpour faire l'ajustement et summarypour calculer les valeurs de p. Bien que cela soit inefficace, il teste le code réel qui a été utilisé. Nous pouvons encore faire des milliers d'itérations en une seconde:

n.sim <- 1e3
set.seed(17)
data.simulated <- matrix(rnorm(n.sim*length(y.expected), y.expected, sigma), ncol=n.sim)
slope.p.value <- function(e) coef(summary(lm(y.expected + e ~ x)))["x", "Pr(>|t|)"]
p.values <- apply(data.simulated, 2, slope.p.value)

Des valeurs de p correctement calculées agiront comme des nombres aléatoires uniformes compris entre et $0$ $1$ lorsque l'hypothèse nulle est vraie. Un histogramme de ces valeurs p nous permettra de vérifier cela visuellement - est-ce qu'il semble à peu près horizontal - et un test d'uniformité khi permettra une évaluation plus formelle. Voici l'histogramme:

h <- hist(p.values, breaks=seq(0, 1, length.out=20))

et, pour ceux qui pourraient imaginer que ce n'est pas suffisamment uniforme, voici le test du chi carré:

chisq.test(h$counts)

X au carré = 13,042, df = 18, valeur p = 0,7891

La grande valeur de p dans ce test montre que ces résultats sont cohérents avec l'uniformité attendue. En d'autres termes, lmc'est correct.

D'où viennent donc les différences de valeurs de p? Vérifions les formules probables qui pourraient être invoquées pour calculer une valeur de p. Dans tous les cas, la statistique de test sera

| t | = | \frac{\hat{β} - 0}{se (\hat{β})} |,

$|t| = \left| \frac{\hat\beta - 0}{\operatorname{se}(\hat \beta)}\right|,$

égal à l'écart entre le coefficient estimé et l'hypothèse (et la valeur correcte) , exprimé en multiple de l'erreur-type de l'estimation du coefficient. Dans la question, ces valeurs sont $\hat \beta$ $\beta = 0$

| t | = | \frac{3.05}{0.87378} | = 3.491

$|t| = \left|\frac{3.05}{0.87378 }\right| = 3.491$

pour l'estimation d'interception et

| t | = | \frac{- 1.38}{0.31906} | = 4.321

$|t| = \left|\frac{-1.38 }{ 0.31906 }\right| = 4.321$

pour l'estimation de la pente. Habituellement, celles-ci seraient comparées à la distribution de Student dont le paramètre des degrés de liberté est (la quantité de données) moins (le nombre de coefficients estimés). Calculons-le pour l'ordonnée à l'origine: $t$ $4$ $2$

pt(-abs(3.05/0.87378), 4-2) * 2

[1] 0.0732

(Ce calcul multiplie la probabilité de Student gauche par car il s'agit d'un test de contre l'alternative bilatérale ) Il est d'accord avec la sortie. $t$ $2$ $H_0:\beta=0$ $H_A:\beta \ne 0$ lm

Un calcul alternatif utiliserait la distribution normale standard pour approximer la distribution de Student . Voyons ce qu'il produit: $t$

pnorm(-abs(3.05/0.87378)) * 2

[1] 0.000482

Effectivement: biglmsuppose que la distribution nulle de la statistique est normale Normal. Quelle est l'erreur? La réexécution de la simulation précédente en utilisant à la place de donne cet histogramme des valeurs de p: $t$ biglmlm

Près de 18% de ces valeurs de p sont inférieures à , un seuil standard de «signification». C'est une énorme erreur. $0.05$

Voici quelques leçons que nous pouvons tirer de cette petite enquête:

N'utilisez pas d'approximations dérivées d'analyses asymptotiques (comme la distribution normale standard) avec de petits ensembles de données.
Connaissez votre logiciel.

whuber
la source

Bonne réponse (+1). Mais vous prenez qui n'est pas vraiment du big data ... Je pense que l'auteur du package a ignoré le petit cas en faveur du cas typique du big data. Il convient toutefois de le signaler dans l'aide pour éviter ces confusions.

n = 4

$n = 4$

n

$n$

epsilone

Pourquoi lm et biglm dans R donnent-ils des valeurs de p différentes pour les mêmes données?

Réponses: