tag:blogger.com,1999:blog-3240271627873788873.post8619163930985097167..comments2024-03-26T06:46:11.752-04:00Comments on Doing Bayesian Data Analysis: p value from likelihood ratio test is not the same as p value from maximum likelihood estimateJohn K. Kruschkehttp://www.blogger.com/profile/17323153789716653784noreply@blogger.comBlogger9125tag:blogger.com,1999:blog-3240271627873788873.post-76877480069601594902014-09-18T15:45:45.286-04:002014-09-18T15:45:45.286-04:00Ah, I see. A solution is emerging! Very interestin...Ah, I see. A solution is emerging! Very interesting.Jake Westfallhttp://jakewestfall.orgnoreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-37781607977250910872014-09-18T09:01:46.970-04:002014-09-18T09:01:46.970-04:00The fact that you had mentioned centering (but not...The fact that you had mentioned centering (but not scaling) x and y, while in my code I'm centering and scaling x and e (but not y) got me wondering about the role of variance. I hadn't thought through it yet, but it sounds like you have. I look forward to seeing the new pictures.Noah Motionhttps://www.blogger.com/profile/00150446498549219747noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-43638606059694185652014-09-18T08:18:38.932-04:002014-09-18T08:18:38.932-04:00The key conceptual difference between the sampling...The key conceptual difference between the sampling distributions (of MLE beta1 vs G2) is not how they treat the intercept, it's how they treat the variance. In some random samples, MLE beta1 can be large but with large MLE sigma, in which case G2 would be relatively small -- these are the blue dots. In other random samples, MLE beta1 can be small but with small MLE sigma, in which case G2 could be relatively large -- these are the red dots. G2 is (in this normal+linear case like the F ratio) a measure of slope vs intercept-only relative to variance in the sample, while MLE beta1 is a 'direct' measure of slope on the absolute scale of the null hypothesis. I'll post some pictures soon...John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-18141881290759171712014-09-18T08:04:41.858-04:002014-09-18T08:04:41.858-04:00I meant to mention this in my first comment, but w...I meant to mention this in my first comment, but when you write that "...it adds another reason never to talk about "the" p value for a set of data, because any data set has many different p values", I feel like it's worth pointing out that anyone who claims that there's a single p value for a data set is forgetting (or doesn't know) that a p value is associated with a <b>test</b> carried out on a data set. A data set without a test doesn't have any p value at all (nor does a test without a data set).<br /><br />Anyway, when I use standardized data (x and e), I get exactly the same p values for the beta coefficient test and the model comparison test. Each time I run 5000 simulations, I get simulated p values that are very close to the values returned by lm() and anova(). I've tried it four or five times with different values for b1 this morning, and the results are very consistent (with respect to the p values matching across tests).<br /><br />Here's my code:<br /><br />N = 100<br />x = runif(N)*5 + 3<br />x = (x-mean(x))/sqrt(((N-1)/N)*var(x))<br />b0 = 0<br />b1 = .23<br />e = rnorm(N,mean=0,sd=7)<br />e = (e-mean(e))/sqrt(((N-1)/N)*var(e))<br />y = b0 + b1*x + e<br /><br />data = data.frame(y=y,x=x)<br /><br />BF = matrix(nrow=nsmp,ncol=2)<br />colnames(BF) = c("beta","Fs")<br />fit1 = lm(y~x,data=data)<br />fit2 = lm(y~1,data=data)<br />test = anova(fit2,fit1)<br /><br />b1.f = coef(fit1)[2]<br />Fs.f = test$F[2]<br /><br />nsmp = 5000<br />b1.H0 = 0<br /><br />for(smpi in 1:nsmp){<br /> e = rnorm(N,mean=0,sd=7)<br /> e = (e-mean(e))/sqrt(((N-1)/N)*var(e))<br /> y = b0 + b1.H0*x + e<br /> d = data.frame(x=x,y=y)<br /> fit.a = lm(y~x,data=d)<br /> fit.b = lm(y~1,data=d)<br /> test.ab = anova(fit.b,fit.a)<br /> BF[smpi,"beta"] = coef(fit.a)[2]<br /> BF[smpi,"Fs"] = test.ab$F[2]<br />}<br /><br />par(mfrow=c(2,1))<br />blims = c(min(-1.25*b1.f,min(BF[,'beta'])),max(-1.25*b1.f,max(BF[,'beta'])))<br />temp = hist(BF[,1],breaks=seq(blims[1],blims[2],length=50))<br />lines(c(b1.f,b1.f),c(0,max(temp$counts)),lwd=3,col="red3")<br />Flim = max(1.25*Fs.f,max(BF[,"Fs"]))<br />temp = hist(BF[,2],breaks=seq(0,Flim,length=50))<br />lines(c(Fs.f,Fs.f),c(0,max(temp$counts)),lwd=3,col="red3")<br /><br />ff = summary(fit1)$fstatistic<br />c(pf(ff[1],ff[2],ff[3],lower.tail=F),2*sum(BF[,1] > b1.f)/nsmp)<br /><br />c(test[["Pr(>F)"]][2],sum(BF[,2] > Fs.f)/nsmp)<br />Noah Motionhttps://www.blogger.com/profile/00150446498549219747noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-41747106282788133372014-09-18T00:14:54.049-04:002014-09-18T00:14:54.049-04:00The joint distribution plot seems to make it prett...The joint distribution plot seems to make it pretty clear that the majority of the discrepancy is due to the one-tailed vs. two-tailed issue that you mentioned in the last post. As for the non-trivial discrepancy that reminds after that, it seems like Noah's suggestion should do it. But you say that when you do this, contrary to Noah, you "get the same results as in the post." But does this include also correcting the tail issue? Between the two measures it seems like we should be able to get them to match pretty much exactly. Having the code could be helpful here.Jake Westfallhttp://jakewestfall.orgnoreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-30778755635726590452014-09-17T23:31:12.698-04:002014-09-17T23:31:12.698-04:00This comment has been removed by the author.John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-73710754851554365522014-09-17T23:31:05.129-04:002014-09-17T23:31:05.129-04:00Noah: When I rerun the simulation with mean-center...Noah: When I rerun the simulation with mean-centered actual data (both x and y), and mean-centered simulated sample data --so that the estimated intercept in actual and sampled data is always zero-- I get the same results as in the post.John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-65541027259251716272014-09-17T21:43:29.389-04:002014-09-17T21:43:29.389-04:00The test of the slope parameter ignores the interc...The test of the slope parameter ignores the intercept parameter and the fact that the estimate will vary depending on whether or not the slope is present in the model.<br /><br />The likelihood ratio test compares the full models, one of which has a slope and an intercept, one of which just has an intercept.<br /><br />I wrote some code to run similar simulations, but first I standardized x and the error before generating y. In this case, the full and restricted models both have estimated intercepts of zero (or close enough, in the ball park of 1e-16), as expected. Both tests also give the same p value (I used the anova() function, so it's a p value for an F statistic, but that should behave pretty much like a G2 statistic).<br /><br />If I plot the joint distribution of b1 and F from the null model, it looks like f(x) = x^2, pretty much - a nice, clean U shape with little to no detectable noise.Noah Motionhttps://www.blogger.com/profile/00150446498549219747noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-38287918517979988402014-09-17T20:43:10.134-04:002014-09-17T20:43:10.134-04:00This recent paper (and its citations) could be use...This recent paper (and its citations) could be useful:<br /><br />http://ba.stat.cmu.edu/journal/forthcoming/smith.pdfAnonymousnoreply@blogger.com