2
$\begingroup$

I'm mainly interested in an estimation of smoothed yearly trend of rate of patients by a disease type. Currently I'm using model with varying coefficients in a form:

    fit <- mgcv::gam(no_patients ~ 
                              s(year, by=disease_type) + 
                              s(clinic_ID, year, bs="re")
                              sex + 
                              age + 
                              disease_type + 
                              offset(population), 
                              family="poisson", 
                              method="REML")

Is there a possibility to get at the same time from the model one (collective) smoothed yearly trend for all disease types?

$\endgroup$

1 Answer 1

2
$\begingroup$

I think you would need to use something like

gam(no_patients ~  s(year) +
                   s(year, by = disease_type, m = 1) + 
                   s(clinic_ID, bs = "re") +
                   sex + 
                   age + 
                   disease_type + 
                   offset(population), 
    family = "poisson", 
    method = "REML")

where we have what we called a "global" smoother s(year) and then specify the factor by smooth to have a 1 derivative penalty (m = 1) in order to allow better identifiability of these repeated smooths of year. This s(year) term is in some sense the smooth of time averaged over all the data (for some definition of "all" - it is still conditional upon all the other terms in the model).

I'm not sure what you were hoping to induce by the s(clinic_ID, year, bs = "re") term? This is going to give you a random linear slope of year for each clinic but not allow each clinic to have it's own intercept. I have removed this term and left only the random intercept per clinic.

If you have more than a few diseases, and you expect the trends per disease are in some way similarly smooth (same amount of wiggliness) then you could use an alternative factor smooth interaction, the fs basis:

gam(no_patients ~  s(year) +
                   s(year, disease_type, bs = "fs") + 
                   s(clinic_ID, bs = "re") +
                   sex + 
                   age +
                   offset(population), 
    family = "poisson", 
    method = "REML")

where we remove m = 1 to allow a second derivative penalty as the fs basis is fully penalised. We also remove the disease_type parametric term because the fs basis contains a random intercept per disease_type.

For more on these models, see our Open Access paper on HGAMs in PeerJ (Pedersen et al 2019).

Pedersen, E.J., Miller, D.L., Simpson, G.L., Ross, N. (2019) Hierarchical generalized additive models: an introduction with mgcv. PeerJ 7:e6876.

$\endgroup$
3
  • $\begingroup$ Thank you. I wanted to include random slope for my clinics as I assume that they have different dynamics not only intercepts. I always thought that random slope is an extension of random intercept. I have only couple diseases so no factor smooth is needed for me. Unfortunately both of your propositions produce singular concurvity in my model and plots for disease types by year say nothing (it doesn't look reasonable, or maybe I should comibine them with the global smoother?). As i noticed in your paper (which is a nice job btw) you just took an average from the smooths - can't I just do it? $\endgroup$
    – Tom
    Sep 10, 2022 at 9:49
  • 1
    $\begingroup$ Could you add outputs from your trials using my suggestions to your original question [edit them in]? I don't understand what "singular concurvity" means - you are combining terms, but the first (having a singular fit) would be worse than having highly concurved smooth functions. What does "say nothing" mean? In and of itself this doesn't seem to be a problem, especially so if the data don't support separate different smooths per disease. I'm not certain now where we used the average of the smooths; perhaps it was to show what the global smooth was doing? 1/2 $\endgroup$ Sep 12, 2022 at 8:11
  • 1
    $\begingroup$ You can always compute something - it's just math after all :-) The key thing will be to estimate the uncertainty of the thing you compute (which you could do using posterior simulation and taking repeated averages of the posterior draws of the smooths) 2/2 $\endgroup$ Sep 12, 2022 at 8:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.