fix for issue with link not adding intercept for linear models

[mattelisi]

We found that there was an issue with the link() function not returning the correct predictions. For example:

library(rethinking)
data(Howell1)
d<-Howell1
d2<-d[d$age>=18,]
m4.3 <- map(
  alist(
    height ~ dnorm( a + b*weight , sigma ) ,
    a ~ dnorm( 178 , 100) ,
    b ~ dnorm( 0 , 1) ,
    sigma ~ dunif( 0 , 50 )
  ),
  data=d2 )

mu <- link(m4.3, data=d2)
str(mu)
 num [1:1000, 1:352] 45.2 42.3 41 41.7 43.6 ...

Note that link give the predictions that are far from the actual values of heights (136+)

I tracked this down to the use of as.character to get formula of the model

> as.character( m4.3@formula[[1]][[3]][[2]] )
[1] "+"          "a"          "b * weight"

I think this is related to the formula being a language type of object. I read in the help that deparse, which is normally preferable to as.character for language objects. Indeed using deparse :

> deparse( m4.3@formula[[1]][[3]][[2]] )
[1] "a + b * weight"

After loading the updated version (devtools::install_github("mattelisi/rethinking", ref="test-link-issues")) the numbers are correct

> mu <- link(m4.3, data=d2)
> str(mu)
 num [1:1000, 1:352] 157 157 157 157 158 ...

An alternative approach that seems to work is as.character( c( m4.3@formula[[1]][[3]][[2]] )).

[rmcelreath]

The issue here is that the model doesn't have a linear model. The version of m4.3 in the main text includes an explicit linear model, and I believe that works, yes?

[mattelisi]

Ah! Yes indeed, I just tried and it works.

I think the issue is then that map() returns correct results regardless of whether the linear formula is inside or outside of dnorm(), no errors or warnings are returned, but then link() gives different results depending on whether the linear formula is inside or outside the density. The documentation of map does not explicitly warn against putting the linear formula within the density function, so some of our students put it inside and we were unsure about why they were getting correct estimates of parameters but wrong predictions.

Switching to deparse, or adding c() in the part that extract the model formula seems to give correct results regardless of whether the linear formula is inside or outside the density.

[rmcelreath]

Okay, thanks. I will think about whether link() on a model without an explicit linear model should fail loudly. Or I could fix it, using your approach, but I'd need some new tests.

Your note about using c() makes it clear what is happening, I think. So thanks!

[rmcelreath]

I pulled this manually into the experimental branch. 3941b14

I still need to run all the test before I merge into master. And I need some new tests for this specific issue now. But I think it works.