I am trying to fit a poisson GLM to count data.

Where the count $N_r$ for time $t_r$ is assumed to follow $\mathrm{ Poisson} (\lambda (t_r - t_0)^{-\alpha} e^{-\beta t_r})$.

Here $r$ is discrete and we observe $N_r$ and $t_r$ for $r$ from 1 to 30. The parameters that I want to fit are $\lambda, \alpha, \beta, t_0 $. Intuitively I am trying to see in what form the counts $N_r$ decrease as time $t_r$ increases.

In R, I can fit the model as follows in the absence of $t_0$ (this is not what I wanted to do, but what i have been able to do) as follows.
$$ model1 = glm(N_r \sim t_r + \log(t_r), family = poisson()) $$ Since, $\log(E(N_r))$ is a linear combination $\log(\lambda)$, $t_r$, and $\log(t_r)$.

The results of model1 is reasonable.

However, I am having difficulty when trying to include the $t_0$ term in my model equation. I would like to fit a model that would estimate the parameters $\lambda, \alpha, \beta, t_0 $. How can i do this in R? It is confusing for me because it appears that, $$ \log(E(N_r)) = \log(\lambda) - \beta t_r - \alpha \log(t_r - t_0)$$ So $\alpha$ and $t_0$ are not separated. (Perhaps I have to use a different function?) It would be great to get some insight on how to fit this model.


2 Answers 2


Unless you have just the right combination of $\alpha$, $\beta$, and values of $x_i$, it may be impossible to pin down the value of $t_0$. The problem is that these curves can be difficult to distinguish.

To obtain a fit, consider searching over a reasonable set of values of $t_0$. For each value the problem is a standard Poisson Generalized Linear Model, and so can be fit routinely. Extract a measure of the fit, such as deviance (or equivalent measures AIC and BIC). A "plausible" value of $t_0$ is any one whose fit isn't too much greater than the value with smallest deviance. A good cutoff threshold to use for the difference is twice an upper percentile of the chi-squared distribution with 1 degree of freedom. The range of resulting values of $t_0$ is a confidence interval for $t_0$.

By choosing suitable units of measurement for $t_r$ (which I will simply call $t$), we may suppose the sample values of $t$ range between $0$ to $1$. This makes it easier to compare models and gives $t_0$ an inherent meaning. For instance, $t_0\approx -1$ is fairly large because it's comparable to the full range of the $t$. We will limit our search to negative values of $t_0$, because the model makes no sense when $t_0$ equals or exceeds the smallest $t$ in the dataset.

Here is an example generated from the model with $\log(\lambda)=0$, $\alpha=2$, $\beta=1/3$, and $t.0 = -1/10$, with $30$ $t$ values equally spaced throughout its range, as described in the question.

Figure 1: data and true model

I searched a range of $t.0$ values geometrically spaced from $-0.01$ to $-3$. Here is the plot of how the deviance varies:

![Figure 2: Deviance vs. t_0

Points within half the 95th percentile of $\chi^2(1)$ (equal to $1.92$) are denoted "significant." The corresponding values of $t_0$ range from $-0.022$ to $-0.37$. Such a narrow range of estimates is rare: in many cases the deviance will scarcely change, no matter what value is given to $t_0$.

The next figure helps show why estimates of $t_0$ can range so widely. For each one of the "significant" values of $t_0$ in the search, I plot the curve corresponding to the GLM fit of $\alpha$, $\beta$, and $\log(\lambda)$. You can see the curves hardly vary:

Figure 3: Plausible fits

They vary the most at low values of $t$. This is valuable information. It indicates that if you had a choice concerning which values of $t$ to measure, you should concentrate many of them at small values. As an example, I created a dataset in which $t$ was sampled at $6$ geometrically spaced values from $0.007$ to $1$. Each one of these values was measured $5$ times independently, producing another dataset of $30$ values. Here are the plausible fits along with the data:

![Figure 4: Plausible fits, alternative data

There is almost no visible difference among them. The range of plausible values of $t_0$ is noticeably narrower, demonstrating the improvement in estimation that is afforded by a better choice of values of $t$ to sample.

To obtain confidence intervals for $\alpha$, $\beta$, and $\log(\lambda)$, extract the confidence intervals from the GLM model fit with $t_0$ set to the one with smallest deviance.

The fits manage to remain the same by rebalancing the effects of the other parameters as $t_0$ is varied. Using this first dataset again, here is how the estimates vary with $t_0$ (as $t_0$ ranges through its $95\%$ confidence interval--the "significant" values shown in the second figure).

![Plots of coefficient estimates

(The true value of $t_0$ is marked with a vertical line.)

Even if you restrict the models to non-negative values of $\alpha$ and $\beta$, there's still quite a range of $(\log(\lambda),\alpha,\beta,t_0)$ parameter combinations that produce almost the same fit.

For those interested in pursuing this further, here is the R code that produced the figures.

# Specify the data-generation model.
lambda.log <- 0
t.0 <- -0.1
alpha <- 2
beta <- 1/3
size <- 0.95 # Percentile of chi-square(1) for deviance cutoff
# n <- 6
# x <- rep(exp(seq(-5, 0, length.out=n)), 5)
x <- seq(0, 1, length.out=30)
f <- function(x, theta) {
  exp(theta["lambda.log"] - log(x - theta["t.0"]) * theta["alpha"] - theta["beta"] * x)
theta <- c(lambda.log=lambda.log, alpha=alpha, beta=beta, t.0=t.0)
# Generate data.
X <- data.frame(x=x, y.hat=f(x, theta))
X$y <- with(X, rpois(length(x), y.hat))
# Generate a closely-spaced set of regressors for the purpose of plotting
# fits, etc.
x.0 <- seq(min(x), max(x), length.out=201)
X.0 <- data.frame(x=x.0, y.hat=f(x.0, theta))
# FIGURE 1: Data with the underlying model.
ggplot(X, aes(x,y)) + geom_path(aes(x,y.hat), X.0, size=1) + 
  geom_point(alpha=1/2, shape=21, fill="Black") + #geom_smooth() + 
  xlab("t") + 
  ggtitle("Data with True Curve")
# Search a range of values of t.0.
u <- -exp(seq(-5, 1, length.out=31))
fits <- lapply(u, function(u) glm(y ~ I(log(x-u)) + x, X, family="poisson"))
# Identify which values of t.0 are closest to achieving the best deviance.
deviance <- sapply(fits, function(x) x$deviance)
i <- which(deviance <= min(deviance) + qchisq(size, 1)/2)
# FIGURE 2: Deviance vs. t.0
D <- data.frame(t.0=u, deviance=deviance, 
                Significant=deviance <=  min(deviance) + qchisq(size, 1)/2)
ggplot(D, aes(x=t.0, y=deviance)) + 
  geom_line(size=1) + 
  geom_point(aes(color=Significant, shape=Significant), size=2) + 
  ggtitle("Deviance vs. t.0")
# FIGURE 3: Display plausible fits.
m <- length(i)
X.1 <- data.frame(t.0=rep(u[i], each=length(x.0)),
                  x = rep(x.0, m), 
                  y.hat=c(sapply(i, function(j) {
                    b <- c(lambda.log=1, alpha=-1, beta=-1) * fits[[j]]$coef
                    f(x.0, c(b, t.0=u[j]))})))
X.1$t.0 <- ordered(signif(X.1$t.0, 1))
ggplot(X.1, aes(x,y.hat)) + 
  geom_line(aes(color=t.0, group=t.0), size=1, alpha=1/2) + 
  geom_point(aes(x,y), X, alpha=1/2, shape=21, fill="Black") +
  xlab("t") + ylab("y") + 
  ggtitle("Plausible Fits For Varying t.0")
# FIGURE 4: Plot coefficients vs. t.0.
a <- sapply(fits[i], coef)
A <- as.data.table(t(a))
names(A) <- c("lambda.log", "alpha", "beta")
A$t.0 <- u[i]
A <- rbind(A[, .(t.0, Estimate=lambda.log, Coefficient="log(lambda)")],
           A[, .(t.0, Estimate=-alpha, Coefficient="alpha")],
           A[, .(t.0, Estimate=-beta, Coefficient="beta")])
ggplot(A, aes(t.0, Estimate, group=Coefficient)) + 
  geom_vline(xintercept=t.0) +
  geom_line(aes(color=Coefficient), size=1) + 
  ggtitle("Coefficient Estimates")
  • $\begingroup$ it is possible to have R source code of this analysis? Many thanks. $\endgroup$
    – Maximilian
    Commented Nov 6, 2017 at 11:28
  • 1
    $\begingroup$ @Maximilian I have now posted it at the end. $\endgroup$
    – whuber
    Commented Nov 6, 2017 at 14:06
  • $\begingroup$ great code and explanation! $\endgroup$
    – Brian
    Commented Nov 21, 2017 at 22:47

Even the model code you gave does not do what you seem to say it does. It assumes $$\log \text{E} N_r = \beta_0 + \beta_1 t_r + \beta_2 \log t_r$$ and estimates $\beta_0$, $\beta_1$ and $\beta_2$, which I believe is not what you had in mind.

You may need to use some software that lets you specify non-linear transformations of parameters or lets you specify a completely user defined likelihood. For the latter there are a lot of options: e.g. PROC NLMIXED in SAS, the rstan R package with it's maximum likelihood capabilities (or perhaps if your parameters are not that well identified you could go Bayesian) and probably a lot of other packages.

  • $\begingroup$ Tks Bjorn, the equation you wrote is what i have in mind (for the simple case with out t zero). Except that there would be a $\beta_0$ . As you said entering the $t_0$ would make the linear model not work. I will try look up the rstan package you mentioned. $\endgroup$
    – Brian
    Commented Nov 5, 2017 at 17:19
  • $\begingroup$ Ok. I'll edit in the $\beta_0$ for clarity. $\endgroup$
    – Björn
    Commented Nov 5, 2017 at 17:43

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.