# How to fit a poisson glm in R when the parameters are not separated

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.

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.

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

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:

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:

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).

(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.

library(ggplot2)
#
# 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.
#
set.seed(17)
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. 3 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. # library(data.table) 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")

• it is possible to have R source code of this analysis? Many thanks. – Maximilian Nov 6 '17 at 11:28
• @Maximilian I have now posted it at the end. – whuber Nov 6 '17 at 14:06
• great code and explanation! – Brian Nov 21 '17 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.

• 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. – Brian Nov 5 '17 at 17:19
• Ok. I'll edit in the $\beta_0$ for clarity. – Björn Nov 5 '17 at 17:43