avehtari / ROS-Examples

Regression and other stories R examples
https://avehtari.github.io/ROS-Examples/
325 stars 256 forks source link

Possible errors in code p. 116 #60

Closed nathanihoff closed 4 years ago

nathanihoff commented 4 years ago

On p. 116, since new is a data frame, I had to make a couple of changes for the code to run for me in the two "by hand" examples.

Just before the second header, I believe the code should be: y_linpred <- a + b*as.numeric(new)

Then in the last line of the second code chunk after the second header: y_pred <- a + b*as.numeric(new) + rnorm(n_sims, 0, sigma

avehtari commented 4 years ago

Thanks for reaching out. These "by hand" examples are part of the notebook https://avehtari.github.io/ROS-Examples/ElectionsEconomy/hibbs.html, they match the book and run for me. Can you try running that notebook? I'm trying to figure what is different in your system.

nathanihoff commented 4 years ago

Thanks for the response. Yes the notebook runs for me, but the problem remains. Sorry, I should have been more clear about the issue: R doesn't throw an error, but the "by hand" code does not replicate the rstanarm functions. With the code in the book/notebook, y_linpred is a scalar, instead of a vector of predicted values. And y_pred is a vector, but it is not the right one; since as.numeric(a + b*new) creates a scalar, y_pred is just one linear predicted value of y, plus the "noise" from rnorm(n_sims, 0, sigma), rather than incorporating the uncertainty of the fitted regression line as well.

For the y_pred example, here's a set of density plots comparing what's in the book (both the rstanarm function and the "by hand" function as written) and my suggested fix: plot_zoom

And the code, with a formal test:

set.seed(1885)
y_pred_rstanarm <- posterior_predict(M1, newdata = new)
set.seed(1885)
n_sims <- nrow(sims)
sigma <- sims[,3]
y_pred_byhand <- as.numeric(a + b*new) + rnorm(n_sims, 0, sigma)
set.seed(1885)
y_pred_byhand_fixed <- a + b*as.numeric(new) + rnorm(n_sims, 0, sigma)

par(mfrow = c(3,1))
plot(density(y_pred_rstanarm))
plot(density(y_pred_byhand))
plot(density(y_pred_byhand_fixed))

all.equal(as.numeric(y_pred_rstanarm), y_pred_byhand)

[1] "Mean relative difference: 0.01517228"

all.equal(as.numeric(y_pred_rstanarm), y_pred_byhand_fixed)

[1] TRUE

avehtari commented 4 years ago

Wow! And Thanks! My bad that I didn't check the results, and I'm still too novice with R to expect that it could do that! I'll fix this tomorrow (now too late).

nathanihoff commented 4 years ago

Happy to help! The unpredictable ways R combines objects of different classes/lengths is one of its more annoying features... I'm really enjoying the book, by the way.

avehtari commented 4 years ago

Thanks again! Fixed now for future printings, in notebook in git repo, and added to the errata.