tag:blogger.com,1999:blog-3240271627873788873.post8764685661559223522..comments2024-03-26T06:46:11.752-04:00Comments on Doing Bayesian Data Analysis: Posterior predicted distribution for linear regression in JAGSJohn K. Kruschkehttp://www.blogger.com/profile/17323153789716653784noreply@blogger.comBlogger14125tag:blogger.com,1999:blog-3240271627873788873.post-34101196994698873292016-10-22T12:06:36.871-04:002016-10-22T12:06:36.871-04:00A more recent post extends this to multiple linear...A more recent post extends this to multiple linear regression:<br /><a href="http://doingbayesiandataanalysis.blogspot.com/2016/10/posterior-predictive-distribution-for.html" rel="nofollow"><br />http://doingbayesiandataanalysis.blogspot.com/2016/10/posterior-predictive-distribution-for.html</a>John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-23147542706424584062016-09-07T20:58:29.200-04:002016-09-07T20:58:29.200-04:00I think it's a matter of being careful with te...I think it's a matter of being careful with terminology. You can say correctly that "for x=[], the 95%HDI for the posterior predicted value of y goes from y=[] to y=[]". But saying that predicted y is "credibly different from zero because the 95%HDI does not include zero" is an extra conclusion based on a decision rule that I eschew. Essentially that's like using a ROPE around y=0 with zero width, and ROPEs should have non-zero width. (In the frequentist world, on the other hand, when a value falls outside the 95% confidence interval it is rejected.) So, I would say you should just report the 95%HDI on the posterior predicted value and leave it at that, without extending it to a conclusion about rejecting zero or being "credibly different" from zero.John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-91713487482237978792016-09-06T08:39:40.831-04:002016-09-06T08:39:40.831-04:00Dear Professor Kruschke – Dear all,
Is it legitim...Dear Professor Kruschke – Dear all,<br /><br />Is it legitimate to use the posterior predictions to make inferences about different values of the predictor?<br /><br />For instance, in the above example, can we conclude that the weight variable is credibly different from, say, 225 when x = 50 but not when x = 80? Can we ask this question for every value of x in order to make inferences about x?<br /><br />Cheers<br />Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-40517771963614605602016-07-13T09:57:46.718-04:002016-07-13T09:57:46.718-04:00@Sean S, you're missing a parenthesis at the e...@Sean S, you're missing a parenthesis at the end of line 36Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-28289646390519784522016-05-12T23:46:49.096-04:002016-05-12T23:46:49.096-04:00I've tried setting this for the multiple regre...I've tried setting this for the multiple regression version, but the JAGS model won't compile.<br /><br />The model below gives this error:<br /> Error parsing model file:<br /> syntax error on line 37 near "yP"<br /><br />JAGS model syntax:<br />1 | data{<br />2 | ym <- mean(y)<br />3 | ysd <- sd(y)<br />4 | for ( i in 1:Ntotal ) {<br />5 | zy[i] <- ( y[i] - ym ) / ysd<br />6 | }<br />7 | for ( j in 1:Nx ) {<br />8 | xm[j] <- mean(x[,j])<br />9 | xsd[j] <- sd(x[,j])<br />10 | for ( i in 1:Ntotal ) {<br />11 | zx[i,j] <- ( x[i,j] - xm[j] ) / xsd[j]<br />12 | }<br />13 | }<br />14 | Nprobe <- length(xP)<br />15 | for ( j in 1:length(xP) ) {<br />16 | zxP[j] <- ( xP[j] - xm ) / xsd<br />17 | }<br />18 | }<br />19 | model{<br />20 | for ( i in 1:Ntotal ) {<br />21 | zy[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zx[i,1:Nx] ) , 1/zsigma^2 , nu )<br />22 | }<br />23 | # Priors vague on standardized scale:<br />24 | zbeta0 ~ dnorm( 0 , 1/2^2 ) <br />25 | for ( j in 1:Nx ) {<br />26 | zbeta[j] ~ dnorm( 0 , 1/2^2 )<br />27 | }<br />28 | zsigma ~ dunif( 1.0E-5 , 1.0E+1 )<br />29 | nu ~ dexp(1/30.0)<br />30 | # Transform to original scale:<br />31 | beta[1:Nx] <- ( zbeta[1:Nx] / xsd[1:Nx] )*ysd<br />32 | beta0 <- zbeta0*ysd + ym - sum( zbeta[1:Nx] * xm[1:Nx] / xsd[1:Nx] )*ysd<br />33 | sigma <- zsigma*ysd<br />34 | # Predicted y values as xProbe:<br />35 | for ( i in 1:Nprobe ) {<br />36 | zyP[i] ~ dt( zbeta0 + sum( zbeta[1:Nx] * zxP[i,1:Nx] , 1/zsigma^2 , nu )<br />37 | yP[i] <- zyP[i] * ysd + ym<br />38 | }<br />39 | }Sean Shttps://www.blogger.com/profile/03902084229952883298noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-11266087078451355112016-04-29T13:06:29.362-04:002016-04-29T13:06:29.362-04:00OK, thank you.
I was speaking about 9.3 Shrinkage ...OK, thank you.<br />I was speaking about 9.3 Shrinkage in hierarchical models.skanhttps://www.blogger.com/profile/03631114761711061847noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-65001145382656939212016-04-29T12:42:44.507-04:002016-04-29T12:42:44.507-04:00Skan: You don't provide any context for your q...Skan: You don't provide any context for your question, but I suspect you are referring to hierarchical models. When displayed as a diagram, low-level parameters are usually lower in the diagram than high-level parameters. For example, parameters that describe individuals are lower-level than parameters that describe group tendencies. Mathematically, if the likelihood function can be factored like this:<br />p(D|a,b,c) = p(D|a) p(a|b) p(b|c) p(c)<br />then a is a lower-level parameter than b which is a lower-level parameter than c.John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-34031983492162673062016-04-29T12:32:51.326-04:002016-04-29T12:32:51.326-04:00Hello.
In your book, What is the meaning of "...Hello.<br />In your book, What is the meaning of "low-level parameters" ? skanhttps://www.blogger.com/profile/03631114761711061847noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-63426726299127141182015-10-31T15:44:20.657-04:002015-10-31T15:44:20.657-04:00Info about fanplot can be found here:
https://gja...Info about fanplot can be found here:<br /><br /><a href="https://gjabel.wordpress.com/2012/08/13/the-fanplot-package-for-r/" rel="nofollow">https://gjabel.wordpress.com/2012/08/13/the-fanplot-package-for-r/</a>John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-80993077686887768802015-10-29T20:17:39.928-04:002015-10-29T20:17:39.928-04:00Dear John, dear all,
I have been plotting posteri...Dear John, dear all,<br /><br />I have been plotting posterior predictives using the fanplot package. For y values at different levels (not for regression but for a certain psychometric function in my case; but it could be applied here as well, I think) I densely generate predicted samples also for new levels in between the experimental levels. The higher the certainty, the more defined the course of the function gets. For low certainty, it will be very fuzzy. It's maybe not as accurate as plotting the little sideways gauss curves at some data points, but it's quite intuitive. Here is an example. <br /><br />http://homepages.uni-paderborn.de/jeti/ppexample.png<br /><br />I'd be interested in your opinions regarding the usefulness of this visualization.<br /><br />Best regards<br />JanUnknownhttps://www.blogger.com/profile/17325379970586504373noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-47281001646146543792015-10-28T07:26:53.668-04:002015-10-28T07:26:53.668-04:00Sure. After all, the posterior predictive distribu...Sure. After all, the posterior predictive distribution is merely a function of the posterior parameter distribution, just like a difference or parameters (e.g., μ1-μ2) or an effect size (e.g., μ/σ) is merely a function of parameters. <br /><br />Also be careful if you're going to make a claim about "rejecting" a difference of zero. If you phrase it as a "hypothesis test" regarding the hypothesis that the difference is zero, then you are in the realm of a model comparison and should be considering Bayes factors and I'm not sure how that would work for predictive distributions. Instead if you phrase it in terms of intervals of the posterior predictive distribution then you should use some reasonable ROPE around the difference of zero.John K. Kruschkehttps://www.blogger.com/profile/17323153789716653784noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-62454028451579355932015-10-28T07:17:29.787-04:002015-10-28T07:17:29.787-04:00One question.
In your example for x=50, the mean y...One question.<br />In your example for x=50, the mean y=76.8, and for x=80, the mean y=211. <br />Would it be appropriate, at the same time we are doing the linear regression, to compare these two posterior predictive distributions using a contrast? The question would be whether the model predicts a credibly higher y outcome for x=80 compared to x=50.<br />Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-73053121344849453422015-10-12T04:58:31.037-04:002015-10-12T04:58:31.037-04:00This looks to be very helpful! I'm teaching a ...This looks to be very helpful! I'm teaching a course using your book this semester and know that this will save me hours of work trying to figure out how to do this. <br /><br />Cheers!<br />ToddAnonymoushttps://www.blogger.com/profile/08903264508857624723noreply@blogger.comtag:blogger.com,1999:blog-3240271627873788873.post-12997278954095110562015-10-11T11:36:39.542-04:002015-10-11T11:36:39.542-04:00This is very helpful. Thank you.
I managed to do ...This is very helpful. Thank you.<br /><br />I managed to do this in R using based on <a href="http://www.indiana.edu/~kruschke/DoingBayesianDataAnalysis/Programs/SimpleRobustLinearRegressionJags.R" rel="nofollow">this program</a>. Which looks to be from the first edition of the book. <br /><br />The solution you post here seems better to me for the reason you state, but also because JAGS and R have different functions for sampling from a T distribution. My model in R looks something like this:<br /><br /><i>rt(1,df=mcmcMat[i,"nu"])*mcmcMat[i,"sigma"] + mcmcMat[i,"beta0"] + mcmcMat[i,"beta1"] * testX</i><br /><br />This caused confusion for me for a while but eventually drove home the point that R and JAGS don't use the same functions.<br /><br />Thank you again for posting this.Joe C. Eatonnoreply@blogger.com