How to produce brms equivalent of se.fit = TRUE of predict() on a glm() fit?

Modeling
1

How to produce brms equivalent of se.fit = TRUE of predict() on a glm() fit?

https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/predict

2

predict in typical R syntax does not produce what Bayesians typically call predictions, but instead the estimated linear predictor for each observation. Compare to fitted which produce the inverse-link of the the estimated linear predictor for each observation.

Here is some code where I explain and show the differences with the analogous functions for a brms fit.

library(tidyverse)
library(brms)


tibble(treatment = gl(3,3), 
       outcome = gl(3,1,9), 
       counts = c(18,17,15,20,10,20,25,13,12)) -> 
  dat

glm.D93 <- glm(counts ~ outcome + treatment, 
               family = poisson(), 
               data=dat)

#Compare these two:
predict(glm.D93) # on log-scale
#>        1        2        3        4        5        6        7        8 
#> 3.044522 2.590267 2.751535 3.044522 2.590267 2.751535 3.044522 2.590267 
#>        9 
#> 2.751535

fitted(glm.D93)
#>        1        2        3        4        5        6        7        8 
#> 21.00000 13.33333 15.66667 21.00000 13.33333 15.66667 21.00000 13.33333 
#>        9 
#> 15.66667

#>predict(glm.D93, se.fit = TRUE)
#>$fit
#>       1        2        3        4        5        6        7        8        9 
#>3.044522 2.590267 2.751535 3.044522 2.590267 2.751535 3.044522 2.590267 2.751535 
#>
#>$se.fit
#>        1         2         3         4         5         6         7         8         9 
#> 0.1708987 0.1957890 0.1860374 0.1708987 0.1957890 0.1860374 0.1708987 0.1957890 0.1860374 
#>
#> $residual.scale
#> [1] 1

brm.D93 <- brm(counts ~ outcome + treatment, 
               family = poisson(), 
               data=dat)
#> Compiling Stan program...
# # removed due to space

fitted(brm.D93, scale="linear") # analogous to predict.glm
#>       Estimate Est.Error     Q2.5    Q97.5
#>  [1,] 3.028829 0.1688073 2.691985 3.360506
#>  [2,] 2.567569 0.1973914 2.177968 2.944656
#>  [3,] 2.736935 0.1918133 2.350793 3.099940
#>  [4,] 3.023705 0.1729421 2.672905 3.347251
#>  [5,] 2.562446 0.1969014 2.162140 2.930231
#>  [6,] 2.731811 0.1920164 2.343455 3.100480
#>  [7,] 3.026297 0.1748992 2.667494 3.352458
#>  [8,] 2.565038 0.1998677 2.165055 2.936777
#>  [9,] 2.734403 0.1861132 2.362159 3.084970

fitted(brm.D93) # analogous to fitted.glm
#>       Estimate Est.Error      Q2.5    Q97.5
#>  [1,] 20.96866  3.549395 14.760948 28.80377
#>  [2,] 13.28983  2.638798  8.828347 19.00412
#>  [3,] 15.72498  3.026292 10.493892 22.19663
#>  [4,] 20.87380  3.579830 14.481984 28.42448
#>  [5,] 13.21845  2.585492  8.689714 18.73196
#>  [6,] 15.64423  2.998461 10.417169 22.20860
#>  [7,] 20.93522  3.633932 14.403831 28.57287
#>  [8,] 13.26040  2.630922  8.715084 18.85498
#>  [9,] 15.66754  2.912125 10.613846 21.86680

predict(brm.D93) # analogous to a summary of simulate
#>       Estimate Est.Error Q2.5  Q97.5
#>  [1,] 20.90250  5.816180   11 33.025
#>  [2,] 13.22275  4.388673    6 23.000
#>  [3,] 15.67875  4.886872    7 26.000
#>  [4,] 20.84800  5.817633   11 34.000
#>  [5,] 13.27250  4.457153    6 23.000
#>  [6,] 15.63975  4.976662    7 27.000
#>  [7,] 20.90975  5.850099   11 33.000
#>  [8,] 13.32000  4.535222    5 23.000
#>  [9,] 15.60825  4.991243    7 26.000

# like this:
simulate(glm.D93, 
         nsim = 1000) %>% 
  apply(MARGIN = 1, 
        function(x){return(c("mean" = mean(x), 
                             "standard.error"=sd(x)))}) %>% 
  t()
#>     mean standard.error
#> 1 20.929       4.575248
#> 2 13.380       3.852589
#> 3 15.643       3.907021
#> 4 21.021       4.499645
#> 5 13.204       3.851392
#> 6 15.630       4.009636
#> 7 20.941       4.619617
#> 8 13.188       3.675073
#> 9 15.682       3.868312

Biggest reason why the standard error for the glm simulation is smaller than analogous brms prediction is that it doesn’t take into account the uncertainty of the parameters.

1 Like