marginal_effects of a non-linear model with beta distribution

[ImanYZ]

I have a non-linear model:

prior1 <- prior(normal(0, 1), nlpar = "skillCoef") +
          prior(normal(0, 1), nlpar = "earlinessCoef") +
          prior(normal(0, 1), nlpar = "PracticeHoursCoef") +
          prior(normal(0, 1), nlpar = "assignmentsHoursCoef") +
          prior(normal(0, 1), nlpar = "lectureprepHoursCoef")
midterm_model <- brm(bf(midterm_0_1 ~ skillCoef*skill + earlinessCoef*earliness_Log + 
PracticeHoursCoef*Practice_Hours + assignmentsHoursCoef*assignments_Hours + 
lectureprepHoursCoef*lecture_prep_Hours, skillCoef + earlinessCoef + PracticeHoursCoef + 
assignmentsHoursCoef + lectureprepHoursCoef ~ 1, nl = TRUE), data = midterm_data, 
prior = prior1, family = Beta(link = "logit"))

When I run plot(marginal_effects(midterm_model), points = TRUE) I get the following error message: Error in terms.default(model) : no terms component nor attribute

[paul-buerkner]

Thanks for reporting this problem. Do you have a reproducible example for me?

[ImanYZ]

Thank you so much for your fast response. Please give me some time to generate an anonymized dataset for you.

[paul-buerkner]

If you can just give me a simulated data set with the variable names you used, we should be good to go.

[ImanYZ]

Sorry for waiting. Here is the anonymized dataset: https://drive.google.com/file/d/1ctnwTDOBII6ZYksMiyVszGdhRXkv_az5/view?usp=sharing

The dependent (outcome) variable is "midterm", which represents grades out of 100. However, because of its specific distribution, I divided all the values by 100 and generated "midterm_0_1" to be able to use family = Beta(link = "logit").

Also, I have two more questions that I'd highly appreciate it if you help me with: 1- How can I use marginal_effects with linear brm? 2- I found this package that plots all the central posterior uncertainty intervals in the same diagram. It's very handy, but it's just plotting the coefficients, so when using family = Beta(link = "logit") those will be on the log odds scale. Is there any other package developed for brms to provide a similar visualization but back-translate onto a mean on the outcome scale?

[paul-buerkner]

While I am working on that fix I have one question for you:

Why are you specfying your model using a non-linear formula although this is clearly just a (generalized) linear model (with beta distribution). You seem to make your life a bit complicated.

[paul-buerkner]

hmm the code works for me. Have you updated brms to the latest github version (2.2.0)?

[ImanYZ]

Let me update and try it again. In terms of your other question, my original model was

midterm_model <- brm(midterm_0_1 ~ skill + earliness_Log + Practice_Hours + assignments_Hours + lecture_prep_Hours, data = midterm_data,
   family = Beta(link = "logit"))

However, I was not able to use marginal_effects(midterm_model) with it. I was searching for non-linear modeling of the data that I came accross marginal_effects.

[paul-buerkner]

marginal_effects should of course work with linear models as well. What do you mean exactly with "I was not able to use..."?

[ImanYZ]

When I used it with the linear model, it returned the same error message: "Error in terms.default(model) : no terms component nor attribute" However, it worked with the non-linear model using a normal family, but after changing the family to family = Beta(link = "logit") I got the error message. I think the issue is not with the linear or non-linear model but with the family. I'm still running the model after updating the package.

[ImanYZ]

After updating the package, when I use family = Beta(link = "logit"), I still get the error: "Error in terms.default(model) : no terms component nor attribute"

[paul-buerkner]

Hmm. Here is the code I used together the data set you send me, which works for me. Can you confirm that this code results in an error for you?

midterm_data <- read.csv("PracticeTool_Data.csv")
prior1 <- prior(normal(0, 1), nlpar = "skillCoef") +
  prior(normal(0, 1), nlpar = "earlinessCoef") +
  prior(normal(0, 1), nlpar = "PracticeHoursCoef") +
  prior(normal(0, 1), nlpar = "assignmentsHoursCoef") +
  prior(normal(0, 1), nlpar = "lectureprepHoursCoef")
midterm_model <- brm(
  bf(midterm_0_1 ~ skillCoef*skill + earlinessCoef*earliness_Log + 
     PracticeHoursCoef*Practice_Hours + assignmentsHoursCoef*assignments_Hours + 
     lectureprepHoursCoef*lecture_prep_Hours, skillCoef + earlinessCoef + PracticeHoursCoef + 
     assignmentsHoursCoef + lectureprepHoursCoef ~ 1, nl = TRUE
  ), 
  data = midterm_data, 
  prior = prior1, family = Beta(link = "logit"),
  chains = 2, cores = 2, iter = 1000
)
plot(marginal_effects(midterm_model), points = TRUE)

[ImanYZ]

I'm testing it but I think the only difference between this one and what I tested before is chains = 2, cores = 2, iter = 1000.

[ImanYZ]

I got the same error message "Error in terms.default(model) : no terms component nor attribute".

[ImanYZ]

By the way, sessionInfo() shows that I have the up-to-date version of "brms_2.2.0".

[paul-buerkner]

Well, that sucks, since I can't fix it when it works for me. Can you post the whole sessionInfo(). Also could you try running my code above in a fresh and clean R session?

[ImanYZ]

R version 3.4.4 (2018-03-15) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Sierra 10.12.6

Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] margins_0.3.0 sjmisc_2.7.0 sjlabelled_1.0.8 sjPlot_2.4.1
[5] betareg_3.1-0 Hmisc_4.0-3 Formula_1.2-2 survival_2.41-3
[9] lattice_0.20-35 brms_2.2.0 Rcpp_0.12.16 tidybayes_0.12.0.000002 [13] modelr_0.1.1 ggplot2_2.2.1 shiny_1.0.5 tidyr_0.8.0
[17] dplyr_0.7.4 statsr_0.0.1 devtools_1.13.5

loaded via a namespace (and not attached): [1] backports_1.1.1 blme_1.0-4 plyr_1.8.4 igraph_1.1.2
[5] lazyeval_0.2.1 TMB_1.7.12 splines_3.4.4 svUnit_0.7-12
[9] crosstalk_1.0.0 TH.data_1.0-8 rstantools_1.4.0 inline_0.3.14
[13] digest_0.6.12 htmltools_0.3.6 rethinking_1.59 rsconnect_0.8.8
[17] magrittr_1.5 checkmate_1.8.5 memoise_1.1.0 cluster_2.0.6
[21] matrixStats_0.52.2 xts_0.10-0 sandwich_2.4-0 colorspace_1.3-2
[25] haven_1.1.0 crayon_1.3.4 lme4_1.1-14 bindr_0.1
[29] zoo_1.8-0 glue_1.2.0 gtable_0.2.0 emmeans_1.1.2
[33] sjstats_0.14.1 rstan_2.17.3 abind_1.4-5 scales_0.5.0
[37] mvtnorm_1.0-7 ggeffects_0.3.1 miniUI_0.1.1 xtable_1.8-2
[41] merTools_0.3.0 htmlTable_1.11.0 ggstance_0.3 foreign_0.8-69
[45] stats4_3.4.4 prediction_0.2.0 StanHeaders_2.17.2 survey_3.32-1
[49] DT_0.4 htmlwidgets_1.0 threejs_0.3.1 arrayhelpers_1.0-20160527 [53] RColorBrewer_1.1-2 acepack_1.4.1 modeltools_0.2-21 pkgconfig_2.0.1
[57] loo_2.0.0 flexmix_2.3-14 nnet_7.3-12 tidyselect_0.2.4
[61] rlang_0.2.0 reshape2_1.4.3 munsell_0.4.3 tools_3.4.4
[65] LaplacesDemon_16.1.0 cli_1.0.0 broom_0.4.4 ggridges_0.4.1
[69] stringr_1.3.0 arm_1.9-3 yaml_2.1.18 knitr_1.20
[73] purrr_0.2.4 bindrcpp_0.2 coin_1.2-2 nlme_3.1-131.1
[77] mime_0.5 rstanarm_2.17.4 compiler_3.4.4 bayesplot_1.5.0
[81] shinythemes_1.1.1 rstudioapi_0.7 tibble_1.4.2 stringi_1.1.7
[85] Brobdingnag_1.2-5 forcats_0.3.0 Matrix_1.2-14 psych_1.7.8
[89] nloptr_1.0.4 markdown_0.8 shinyjs_1.0 effects_4.0-0
[93] stringdist_0.9.4.7 pillar_1.2.1 pwr_1.2-2 lmtest_0.9-35
[97] bridgesampling_0.4-0 estimability_1.3 data.table_1.10.4-3 httpuv_1.3.5
[101] R6_2.2.2 latticeExtra_0.6-28 gridExtra_2.3 codetools_0.2-15
[105] colourpicker_1.0 MASS_7.3-49 gtools_3.5.0 assertthat_0.2.0
[109] withr_2.1.0 shinystan_2.4.0 mnormt_1.5-5 multcomp_1.4-8
[113] parallel_3.4.4 grid_3.4.4 rpart_4.1-13 coda_0.19-1
[117] glmmTMB_0.2.0 minqa_1.2.4 snakecase_0.9.0 carData_3.0-0
[121] base64enc_0.1-3 dygraphs_1.1.1.4

[paul-buerkner]

Thanks! I see quite a few packages being attached. Maybe there is some function masking happening somewhere.

Could you run the model in a fresh R session with just brms loaded via library(brms) and no other objects in the global environment.

[ImanYZ]

Sure.

[ImanYZ]

You're right. The issue was because of some function masking. It works fine in a new session. Thank you so much and sorry for taking your time a lot.

Whenever you have some time, I'd appreciate it if you answer my second question as well: I found this package that plots all the central posterior uncertainty intervals in the same diagram. It's very handy, but it's just plotting the coefficients, so when using family = Beta(link = "logit") those will be on the log odds scale. Is there any other package developed for brms to provide a similar visualization but back-translate onto a mean on the outcome scale? The reason I'm looking for this is that in behavioral economics and experimental psychology, we're usually interested in comparing the effect size of different factors under the study on an outcome, e.g., in the current study that I've shared the dataset with you, we are trying to compare the effect of spending an hour on different learning activities (Practice_Hours, assignments_Hours, lecture_prep_Hours) on the midterm grade. Having the confidence intervals in odds ratio is very difficult to interpret.

[paul-buerkner]

You are welcome. With regard to your second questions, I suggest you take a look at the tidybayes package (https://github.com/mjskay/tidybayes), which works nicely together with brms.

[ImanYZ]

Thank you again.