# Improve estimate of zeros and how to deal with concurvity in GAM NB

I'm trying to model the number of deaths for each age, for 4 regions, across 5 years. I fitted several models, but end up with the following model:

m1 <- gam(deaths ~ ti(region, bs='re') + ti(ageCenter, bs='cr', k=40) +
ti(ageCenter, region, k=c(40,4)) + ageCenter*yearCenter + offset(log(PopMedia)),
family = nb, data = male, method = "REML")


The model is reasonable, all variables are significant. However, the concurvity is not great for region:

         para ti(Regiao) ti(ageCenter) ti(ageCenter,Regiao)
worst       1  0.8335333   1.000000000         0.0635483181
observed    1  0.8335333   0.007150858         0.0003614262
estimate    1  0.8335333   0.029291067         0.0005498903


Nevertheless, I compared the observed number of deaths with the predicted and the estimate was low for one region (the smaller region). I tried to improve the model and I get a closer estimate for the smaller region but another region decreased the estimate number of zeros. This is the new model:

m2 <- gam(Obitosdx ~ ti(Regiao, bs='re') + ti(ageCenter, bs='cr', k=40, by=Regiao) +
ti(ageCenter, Regiao, k=c(40,4)) + ageCenter*yearCenter + offset(log(PopMedia)),
family = nb, data = male, method = "REML")


All variables are still significant, the concurvity remains unchanged as there's still ti(region) in the model. Residuals between the two models are similar.

My first question relates to the models. I'm not sure if it's appropriate to have ti(ageCenter, bs='cr', k=40, by=Regiao) in the same model as ti(ageCenter, Regiao, k=c(40,4)). Wouldn't the ti product interaction also give different smooths for each level of Region?

My second question is the concurvity for ti(region) worrisome? If the variable is significant how to deal with this?

These are the two models, comparing the m1 and m2 for the smallest region and a region intermediate. Do you have any suggestion for improvement? The estimate of the number of zeros is quite low. For the smallest region, the m2 estimates a closer number of zeros, while for the intermediate region, the m1 estimates a closer number of zeros. Any suggestion is most welcome! EDIT

Following @Gavin suggestion, I've recoded Region to a factor variable, which should have been a factor to begin with. However, I still get a large deviation for young ages. Any suggestion? The plot above with the intermediate region corresponds to Alentejo, while the smallest to RAM. You can't fit a smooth to a factor variable and this means that your ti(ageCenter, Regiao, k=c(40,4)), ti(Regiao, bs = "re"), and ti(ageCentre, by = Regiao) are not doing what you might expect them too.

That the former even works at all indicates that Region (or Regiao, which is it? Does it matter?) is not stored as factor variable. That means that the by smooth is actually a linear effect of Region modified by a smooth of ageCentre, and not separate smooths for the regions.

The high concurvity of the terms is most likely due to the above problems.

You are also going to need to be very careful with the ageCentre*yearCentre term as this will include a linear effect of ageCentre and somewhere among the bases of the splines you'll have a linear effect of ageCentre also, somewhere in the null space of one of the splines. This is going to introduce some identifiability issues; or it may, unless mgcv is being very clever to remove those effects (and I don't think it is doing this).

Assuming Region is a factor (it is isn't, recode it to one), then a model with separate smooths of ageCentre for each region plus the effect of ageCentre and yearCentre might be written as

m1 <- gam(deaths ~ ti(yearCenter) +
region +
ti(ageCenter, bs='cr', k=40, by = region) +
ti(ageCenter, yearCenter) +
offset(log(PopMedia)),
family = nb, data = male, method = "REML")


I made the ageCentre by yearCentre thing a smooth, and I doubt you want to estimate a random effect with four levels (but if you must, put the re basis and smooth back on the region term).

• Thank you for your response! So the problem is probably silly, since aparently Region is not a factor. But it should be =/ I thought I only needed ti(ageCenter, bs='cr', k=40, by = region) but I was confused as to why the model improved with ti(Regiao, bs = "re") and ti(ageCenter, Regiao, k=c(40,4)). I may be misunderstanding something but ageCenter*yearCenter should be similar to a non-smoothed ti? Completely wrong? I didn't included a smooth for year because I didn't thought it was relevant but an interaction age year would be, as we will tend to have an older population. – psoares May 30 '18 at 18:01
• Thank you for your answer @Gavin. It was a stupid mistake. However, I still get a large deviation for the young ages, as seen in the residuals figure uploaded. I think an adaptive smooth would work instead of a ti. Do you think that would be appropriate? – psoares Jun 1 '18 at 7:32
• You may have an issue in that the model may not be flexible enough to simultaneously predict low values and some higher counts. This would imply the expected count varies more than is allowed under the estimated NB distribution. This could be because the overdispersion parameter is not constant for all observations but may itself vary with covariates (i.e. you want a model for theta), or you may have zero-inflation, or the NB may not be a good conditional distribution for the response for other reasons. – Gavin Simpson Jun 1 '18 at 15:45
• I would look at a rootogram of the residuals (see a post i wrote on my blog for details on how to do that: fromthebottomoftheheap.net/2016/06/07/rootograms ) – Gavin Simpson Jun 1 '18 at 15:45
• I actually think the GAM with an adaptive smooth works well. This model gam(deaths ~ yearCenter + Region + s(ageCenter, bs='ad', k=50, by=Region) + offset(log(PopMedia)), family = nb, data =female, method = "REML") is able to predict a closer number of zeros than using the ti. I will have a look at your page. Thank you for your help! – psoares Jun 1 '18 at 15:56