I am not sure if I am using the correct terminology, something must be written about the following problem, but I cannot find it by searching.

I am presently analyzing data about the effect of introduction of a new vaccination program. The data is supposedly "ecological data", that is , population data, we do not have data on individual level, so no subject level covariates. So, basically, I have a count time series, of total number of deaths, and number of deaths for reason A. But we do not know how good the surveillance system is, that is, which percentage of cases is actually reported. That percentage is surely changing with time, hopefully it is becoming better. But we do assume, realistically, that the probability of a death being reported do not depend on its cause.

So, I model $$ Y_{1t} \sim \text{Po}(\lambda_t) \\ Y_{2t} \sim \text{Po}(\mu_t) $$ where $Y_{1t}$ is number of deaths for reason A at $t$, while $Y_{2t}$ is number for all other reasons. In reality, those variables are overdispersed, but here I will just assume Poisson (because overdispersion issues are unimportant for my question).

So, I use glms (actually, vglms from R package VGAM) for representing Poisson regressions, with some common terms in the linear predictors $\eta_1, \eta_2$. I introduce a 0/1 dummy for the introduction of the vaccine at time $t_0$, a polynomial of low degree for representing the observed relative decrease of deaths of certain reason relative for other reasons, a common natural spline term for representing the effect of the varying coverage (assumed the same for both processes) and some other terms unimportant for the question. The model is then more explicit: $$ Y_{1t} \sim \text{Po}(\lambda_t=\exp(\beta_{10}+\beta_{11}t+\beta_{12}t^2+\beta_{13}t^3+s(t)+\beta_{14}I(t-t_0\ge 0)) ) \\ Y_{2t} \sim \text{Po}(\mu_t=exp(\beta_{20}+s(t))) $$ where $s(t)$ is a common spline term. Since the two regressions have some common parameters, we use VGAM for the estimation.

My problem is that by introducing too much flexibility (too many df for the natural spline term), that term will eventually absorb the effect of the dummy term! So, how to choose how much flexibility? As we do not really have prior info about the correct nuber of df, it must in some way be from the fit.

Some observations: Since the smooth term is common for both linear predictors, if too little flexibility, will cause underfitting, which can be seen as autocorrelated residuals, and even correlation between the two residuals. So that will give some information.

But, I am only concerned about good estimation (unbiased, small bias, consistency) of the treatment effect. How can I be sure about that, in the presence of many other badly identified parameters, with little prior information? That is, should I underfit or overfit (my intuition says that I should avoid overfitting, so erring in the direction of some underfit) but there must be some papers studying this or some related questions?

Any thoughts, or references?


I found one paper which seems relevant, http://ac.els-cdn.com/S0304414913000811/1-s2.0-S0304414913000811-main.pdf?_tid=5ddf25ce-ac17-11e6-a0a4-00000aacb35d&acdnat=1479312884_12ea86a56c9581553b7607211e21941d but there must be something more specific? Also, very low views in this question, is there something I can do to make the question clearer?


As the observation errors acts multiplicatively on the true deaths, including it additively on the linear scale doesn't seem optimal to me. At least in my field, the standard would be a hierarchical model that first models the (latent) true events for reason A as

expectedDeathsA = exp(x1 + x2 + f(vaccine, t) )

trueDeathA ~ poisson(expectedDeathsA)

and then derive from that the two observed responses

reportedDeathsA ~ binomial(trueDeathA, observationError)
reportedTotalDeaths ~ f(trueDeathA, otherDeaths)

where observationError could be again a function of time if needed, and the "~ f" in reportedTotalDeaths could mean any function or distribution you want. It's pretty easy to implement such a model in JAGS or STAN, but you will likely have to go Bayesian (no idea if AD Model Builder would work as well, if yes this could provide MLEs).

A word of warning: the experience is that it is very hard to estimate the observation accuracy if you don't have a few cases where you know trueDeathA, or a strong prior (but then you don't really estimate the accuracy). Even with repeated observations (which you don't have if the vaccine effectiveness changes in time), small model errors / changes can have huge effects on the estimates. The idea with including the reportedTotalDeaths is helpful to avoid this problem, but I'm not sure how much benefits it generates, as you have to consider that otherDeaths could change over time as well.


  • A standard reference for the model without the totalDeaths would be Royle, J. A. (2004) N-Mixture Models for Estimating Population Size from Spatially Replicated Counts. Biometrics, 60, 108-115.

  • An excellent textbook with code on these models is: Kéry, M. & Royle, J. A. (2015) Applied Hierarchical Modeling in Ecology: Analysis of distribution, abundance and species richness in R and BUGS: Volume 1: Prelude and Static Models. Academic Press

  • The following reference discusses data needs for reliably estimating observation errors Guillera-Arroita, G. (2016) Modelling of species distributions, range dynamics and communities under imperfect detection: advances, challenges and opportunities. Ecography.

These are all for species abundances, but in the end it's all count data I guess.

  • $\begingroup$ Thanks, do you have references to some published papers, maybe from your field? $\endgroup$ – kjetil b halvorsen Nov 17 '16 at 13:24
  • $\begingroup$ added a few refereces $\endgroup$ – Florian Hartig Nov 17 '16 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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