Posterior Predictions from Multiple Equations

[rmcelreath]

Posterior Predictions from Multiple Equations

[rmcelreath]

The function link() just substitutes values into any linear models in the model definition. You provided a value for x1, but not for y1. It isn't going to propagate the prediction yhat1 into the linear model for yhat2, because the model isn't written that way.

But you can just put a dummy linear model in the model definition to get what you are after.

nf <- alist(
  y1 ~ dnorm(yhat1, sigma1),
  y2 ~ dnorm(yhat2, sigma2),

  yhat2pred <- b2*yhat1,
  yhat1 <- b1*x1,
  yhat2 <- b2*y1,

  b1 ~ dnorm( 0 , 10 ),
  b2 ~ dnorm( 0 , 10 ),

  sigma1 ~ dcauchy( 0 , 1 ),
  sigma2 ~ dcauchy( 0 , 1 )

)

n.fit <- map2stan( 
  nf , 
  data=ndf , 
  start=list(b1=1, b2=1, sigma1=10, sigma2=10)
)

a <- link(n.fit, data=list(x1=3))
sapply(a,mean)

The output should be:

yhat2pred     yhat1     yhat2 
13.027757  6.586785  1.977875 

I think the yhat2 value is coming from default values for y1 that are used, since you didn't pass in a value to link.