{"action":"created","issue":{"url":"https://api.github.com/repos/rmcelreath/rethinking/issues/15","labels_url":"https://api.github.com/repos/rmcelreath/rethinking/issues/15/labels{/name}","comments_url":"https://api.github.com/repos/rmcelreath/rethinking/issues/15/comments","events_url":"https://api.github.com/repos/rmcelreath/rethinking/issues/15/events","html_url":"https://github.com/rmcelreath/rethinking/issues/15","id":125241429,"number":15,"title":"Posterior Predictions from Multiple Equations","user":{"login":"jebyrnes","id":757410,"avatar_url":"https://avatars.githubusercontent.com/u/757410?v=3","gravatar_id":"","url":"https://api.github.com/users/jebyrnes","html_url":"https://github.com/jebyrnes","followers_url":"https://api.github.com/users/jebyrnes/followers","following_url":"https://api.github.com/users/jebyrnes/following{/other_user}","gists_url":"https://api.github.com/users/jebyrnes/gists{/gist_id}","starred_url":"https://api.github.com/users/jebyrnes/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/jebyrnes/subscriptions","organizations_url":"https://api.github.com/users/jebyrnes/orgs","repos_url":"https://api.github.com/users/jebyrnes/repos","events_url":"https://api.github.com/users/jebyrnes/events{/privacy}","received_events_url":"https://api.github.com/users/jebyrnes/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":0,"created_at":"2016-01-06T19:12:42Z","updated_at":"2016-01-07T08:50:29Z","closed_at":null,"body":"So, I'm working with rethinking in a Structural Equation Modeling framework. It works beautifully, but, I'm having some issues in generating posterior predictions. Namely, effects of linked relationships don't seem to translate.  Consider the following model\r\n\r\nx1 -> y1 -> y2\r\n\r\nwhere each coefficient is a slope of 2.\r\n\r\nI can create a nice fake data frame\r\n\r\n    set.seed(31415)\r\n    ndf <- data.frame(x1=runif(100,0,100))\r\n    ndf$y1 <- rnorm(100, ndf$x1*2, 40)\r\n    ndf$y2 <- rnorm(100, ndf$y1*2, 80)\r\n\r\n\r\nAnd then a model that fits nicely using map2stan\r\n\r\n\r\n    nf <- alist(\r\n      y1 ~ dnorm(yhat1, sigma1),\r\n      y2 ~ dnorm(yhat2, sigma2),\r\n  \r\n      yhat1 <-  b1*x1,\r\n      yhat2 <- b2*y1,\r\n  \r\n      b1 ~ dnorm( 0 , 1000 ),\r\n      b2 ~ dnorm( 0 , 1000 ),\r\n  \r\n      sigma1 ~ dcauchy( 0 , 1 ),\r\n      sigma2 ~ dcauchy( 0 , 1 )\r\n      \r\n    )\r\n\r\n    n.fit <- map2stan( \r\n      nf , \r\n      data=ndf , \r\n      start=list(a1=1, a2=1, b1=1, b2=1, sigma1=10, sigma2=10)\r\n    )\r\n\r\n\r\nhere\r\n\r\n    > coef(n.fit)\r\n            b1         b2     sigma1     sigma2 \r\n      1.930119   1.965780  37.934121  86.538586 \r\n\r\n\r\nGreat!  But, when I try and do posterior prediction with x starting as 3, which should produce roughly 12 for y2...\r\n\r\n\r\n    > link(n.fit, data=list(x1=3)) -> a\r\n    [ 1000 / 1000 ] 200 / 1000 ]\r\n \r\n    > sapply(a, mean)\r\n       yhat1    yhat2 \r\n    5.790357 1.965780 \r\n`\r\n\r\nHuh....why is yhat incorrect here?"},"comment":{"url":"https://api.github.com/repos/rmcelreath/rethinking/issues/comments/169600067","html_url":"https://github.com/rmcelreath/rethinking/issues/15#issuecomment-169600067","issue_url":"https://api.github.com/repos/rmcelreath/rethinking/issues/15","id":169600067,"user":{"login":"rmcelreath","id":3230381,"avatar_url":"https://avatars.githubusercontent.com/u/3230381?v=3","gravatar_id":"","url":"https://api.github.com/users/rmcelreath","html_url":"https://github.com/rmcelreath","followers_url":"https://api.github.com/users/rmcelreath/followers","following_url":"https://api.github.com/users/rmcelreath/following{/other_user}","gists_url":"https://api.github.com/users/rmcelreath/gists{/gist_id}","starred_url":"https://api.github.com/users/rmcelreath/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/rmcelreath/subscriptions","organizations_url":"https://api.github.com/users/rmcelreath/orgs","repos_url":"https://api.github.com/users/rmcelreath/repos","events_url":"https://api.github.com/users/rmcelreath/events{/privacy}","received_events_url":"https://api.github.com/users/rmcelreath/received_events","type":"User","site_admin":false},"created_at":"2016-01-07T08:50:28Z","updated_at":"2016-01-07T08:50:28Z","body":"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.\r\n\r\nBut you can just put a dummy linear model in the model definition to get what you are after.\r\n```\r\nnf <- alist(\r\n  y1 ~ dnorm(yhat1, sigma1),\r\n  y2 ~ dnorm(yhat2, sigma2),\r\n\r\n  yhat2pred <- b2*yhat1,\r\n  yhat1 <- b1*x1,\r\n  yhat2 <- b2*y1,\r\n\r\n  b1 ~ dnorm( 0 , 10 ),\r\n  b2 ~ dnorm( 0 , 10 ),\r\n\r\n  sigma1 ~ dcauchy( 0 , 1 ),\r\n  sigma2 ~ dcauchy( 0 , 1 )\r\n\r\n)\r\n\r\nn.fit <- map2stan( \r\n  nf , \r\n  data=ndf , \r\n  start=list(b1=1, b2=1, sigma1=10, sigma2=10)\r\n)\r\n\r\na <- link(n.fit, data=list(x1=3))\r\nsapply(a,mean)\r\n```\r\nThe output should be:\r\n```\r\nyhat2pred     yhat1     yhat2 \r\n13.027757  6.586785  1.977875 \r\n```\r\nI think the yhat2 value is coming from default values for y1 that are used, since you didn't pass in a value to link."}}