# Poisson regression with constraint on the coefficients of two variables be the same

The aim of this experiment is to explore the effects of age, period and cohort. Thus, none of them can be thrown.

Therefore, the assumption of no cohort effects greatly simplifies estimations but can lead to model misspecification and is inconsistent with accumulating evidence of cohort changes in a variety of health outcomes and mortality'., 'Age-Period-Cohort Analysis'. Yang Yang (2013) P.64

There are numerous approaches to solving the unidentified problem, the one I choose here is constrained generalized linear model, CGLIM. I also tried to add nonlinear term into the model, i.e. $$age^2$$, but NAs still appear.

I tried to use restriktor R package to fit the count data with constrained Poisson regression.

The sample data is generated by the code below:

library(tidyverse)
library(restriktor)

set.seed(123)
period_test = rep(seq(2000, 2020, 5), 3)
age_test = c(seq(30, 50, 5), seq(50, 70, 5), seq(70, 90, 5))
cohort_test = period_test-age_test
period_test = period_test %>% fastDummies::dummy_cols() %>% modify(as.factor)
names(period_test) = names(period_test) %>% gsub('.data_','period_',.)
death_data = c(rpois(5, 30), rpois(5, 20), rpois(5, 20))

testing = data.frame(period_test, age_test, cohort_test, death_data)


Then I fit the data with glm. The age, period and cohort term are unidentified because of linear dependence. Therefore, I would like to use restriktor package to fit the Poisson regression with
equality constraint.

tt = glm(death_data ~ period_2000 + period_2005 + period_2010 +
period_2015 + period_2020 + as.factor(age_test) +
as.factor(cohort_test), family = poisson, data = testing)

my_constraints = 'period_20001 = period_20051'
restriktor(tt, my_constraints)


The glm result for tt is:

Call:  glm(formula = death_data ~ period_2000 + period_2005 +
period_2010 + period_2015 + period_2020 + as.factor(age_test) +
as.factor(cohort_test),
family = poisson, data = testing)

Coefficients:
(Intercept)                period_20001                period_20051                period_20101                period_20151
3.8729                     -0.3466                     -1.2339                     -0.6949                     -0.9826
period_20201       as.factor(age_test)35       as.factor(age_test)40       as.factor(age_test)45       as.factor(age_test)50
NA                      0.6931                      0.3483                      0.6650                     -0.2094
as.factor(age_test)55       as.factor(age_test)60       as.factor(age_test)65       as.factor(age_test)70       as.factor(age_test)75
0.7320                     -0.1865                      0.7943                     -0.2683                          NA
as.factor(age_test)80       as.factor(age_test)85       as.factor(age_test)90  as.factor(cohort_test)1950  as.factor(cohort_test)1970
NA                          NA                     -1.5703                     -0.4266                          NA

Degrees of Freedom: 14 Total (i.e. Null);  0 Residual
Null Deviance:      46.46
Residual Deviance: -1.91e-14    AIC: 104.4


There are NAs in the coefficients. Then, I use the following constraint:my_constraints = 'period_20001 = period_20051'


The result of restriktor function is :
r
> restriktor(tt, my_constraints)
Error in qr.default(t(ceq.JAC)) :
NA/NaN/Inf in foreign function call (arg 1)


What can I do to fit with equality constrained glm ...?

Thank for Dave's reminder, my purpose for fitting constrained Poisson regression is to plot a Forest plot: something like this:

from wiki

Therefore, I am wondering whether the confidence intervals of period_2000, period_2005 will be the same or not if I use Dave's method.

Update

According to the chapterBayesian Age-Period-Cohort Models written by Ethan Fosse, the Classical APC Regression Model is the method I referred to in the question. In my specific case, I opt to set the period effects for the years 2000 and 2005 as equal, imposing constraints during model estimation. The rationale behind selecting these two periods lies in my assumption that their effects are highly similar. Here’s how the implementation will proceed:

set.seed(123)
period_test = rep(seq(2000, 2020, 5), 3)
age_test = c(seq(30, 50, 5), seq(50, 70, 5), seq(70, 90, 5))
cohort_test = period_test-age_test
period_test = period_test %>% fastDummies::dummy_cols() %>% modify(as.factor)
period_test$$period_2000_2005 = period_test$$period_2000 + period_test$$period_2005 period_test$$period_2000=NULL
period_test\$period_2005=NULL
names(period_test) = names(period_test) %>% gsub('.data_','period_',.)
death_data = c(rpois(5, 30), rpois(5, 20), rpois(5, 20))

testing = data.frame(period_test, age_test, cohort_test, death_data)
tt = glm(death_data ~ period_2000_2005 +
period_2010 + period_2015 + period_2020 +
as.factor(age_test) + as.factor(cohort_test),
family = poisson, data = testing)


The result is as below:

Call:  glm(formula = death_data ~ period_2000_2005 + period_2010 +
period_2015 + period_2020 + as.factor(age_test) +
as.factor(cohort_test), family = poisson, data = testing)

Coefficients:
(Intercept)            period_2000_2005                 period_2010
3.66537                    -0.29808                    -0.52988
period_2015                 period_2020       as.factor(age_test)35
-0.52988                          NA                    -0.14842
as.factor(age_test)40       as.factor(age_test)45       as.factor(age_test)50
0.04256                     0.19671                    -0.48732
as.factor(age_test)55       as.factor(age_test)60       as.factor(age_test)65
-0.37609                    -0.44967                    -0.44967
as.factor(age_test)70       as.factor(age_test)75       as.factor(age_test)80
-0.72824                    -0.18924                          NA
as.factor(age_test)85       as.factor(age_test)90  as.factor(cohort_test)1950
NA                    -0.52988                    -0.04676
as.factor(cohort_test)1970
NA

Degrees of Freedom: 14 Total (i.e. Null);  0 Residual
Null Deviance:      16.68
Residual Deviance: -1.554e-14   AIC: 103.1



There is still NA in the result, which indicates that the model has more parameters than can be estimated from the data. So we may need to impose more constraints on the variables.

You may check the chapter written by Ethan Fosse, or the Summer Institutes in Biostatistics of the University of Washington (by Jon Wakefield): Module 3: Age-Period-Cohort Modeling and Analysis Lecture 1: Preliminaries Lecture 2: Identification and Modeling Lecture 3: Splines, Smoothing, Bayes and INLA Lecture 4: Bayesian Methods

Another question related to this topic: Can we use dummy variables to capture the nonlinearity in Age, Period, and Cohort effects?

Can we use dummy variables to capture the nonlinearity in Age, Period, and Cohort effects?

• I cannot get your code to work. Subtraction is not meaningful for factors, and I get a "variable lengths differ" message. I am running on the latest version of all the packages. Specifically, cohort_test = period_test - age_test returns a matrix of all NA values. You have also constructed your dummies such that many of the coefficients are estimated as NA, as you have a lot more variables than observations. Even leaving out cohort_test, you have 18 variables and 15 observations. Commented Jan 18, 2023 at 16:42
• Thank you for pointing out this! I am not familiar with how to construct sample data... Commented Jan 19, 2023 at 2:12
• You will see an example in my updated answer below. Commented Jan 19, 2023 at 2:34
• You are editing an old post with answers to add new questions. That is not how this site is supposed to work, and it is likely to be ineffective, as few people will note the new question. You are better off to ask a new question (with links back to this q!) Commented Feb 27 at 3:38
• Thank you for your reminder! I have asked a new question. Commented Feb 27 at 3:49

Since the cohort equals the period - the age (cohort_test = period_test-age_test,) there certainly is a linear dependence! There is no value to including all three terms in a linear model, and the proper way to address this is to choose two of them and drop the third. If for some reason you want to parameterize the model using a different subset of two of the three terms, you could either rerun the model with the new terms or, with a little more pain, calculate the coefficients of the two "new" terms from the "old" terms.

Another way around this problem is to define new variables based upon the current ones in such a way that there is a nonlinear relationship between them. One technique that preserves the factor-oriented approach above is to group the variables into chunks, e.g., instead of having a factor for an age of 67 and another for an age of 68, group the ages together after calculating the cohort, which can also be grouped. The grouping may well break the linear relationship:

set.seed(123)

period <- rep(seq(2000,2020,5),15)
age = rep(c(seq(31,50,2),seq(51,70,2),seq(71,90,4)), 3)
cohort <- period - age

# block age, cohort into 10 year chunks
age <- 10*(age %/% 10)
cohort <- 10*(cohort %/% 10)

death_data = c(rpois(25,30),rpois(25,20),rpois(25,15))
testing = data.frame(as.factor(period), as.factor(age), as.factor(cohort), death_data)

tt = glm(death_data~., family = poisson,data = testing)
summary(tt)


with result:

Call:
glm(formula = death_data ~ ., family = poisson, data = testing)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-3.4540  -1.4262  -0.2083   1.2389   3.5188

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)             2.4256     0.8177   2.966  0.00301 **
as.factor.period.2005  -0.1349     0.2174  -0.620  0.53497
as.factor.period.2010  -0.1406     0.2174  -0.647  0.51793
as.factor.period.2015  -0.2433     0.2370  -1.027  0.30463
as.factor.period.2020  -0.3316     0.3802  -0.872  0.38318
as.factor.age.40        0.0930     0.2131   0.436  0.66255
as.factor.age.50        0.2152     0.4016   0.536  0.59197
as.factor.age.60        0.3410     0.5929   0.575  0.56521
as.factor.age.70        0.7385     0.8091   0.913  0.36138
as.factor.age.80        0.7694     0.8997   0.855  0.39243
as.factor.cohort.1930   0.2123     0.2578   0.824  0.41015
as.factor.cohort.1940   0.4692     0.4708   0.997  0.31900
as.factor.cohort.1950   0.6600     0.6445   1.024  0.30581
as.factor.cohort.1960   0.7199     0.8296   0.868  0.38552
as.factor.cohort.1970   0.8685     1.0193   0.852  0.39421
as.factor.cohort.1980   0.9817     1.1910   0.824  0.40977
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 219.04  on 74  degrees of freedom
Residual deviance: 211.14  on 59  degrees of freedom


Some trial and error may be required, as small samples for one or two age - period combinations may cause a small number of the cohort factors to be redundant.

• As I mentioned in the question, I tried to add a constraint to solve the unidentified problem: the linear of the effect of period_2000 = period_2005, which is an approach proposed by Mason and Smith. Therefore, what I need is to fit the model with constrained Poisson regression, in order to show the effect of age, period and cohort separately. >~< Commented Jan 18, 2023 at 5:01
• You might want to try: "Nonlinear parametric (algebraic) transformation approach Define a nonlinear parametric function of one of the age, period, or cohort variables so that its relationship to others is nonlinear." Commented Jan 19, 2023 at 1:56
• You only have 15 observations and a lot more variables than that. See my comment in response to your original post... you need to redesign your experiment, maybe make the ages and periods fall into fewer groups and get more data. Commented Jan 19, 2023 at 2:06
• I've expanded my answer with an example above! Commented Jan 19, 2023 at 2:28
• 1. Your dataset still suffers from only having 15 observations. 2. You may have two different ages that fall into the same group, e.g., age 31 and 38, but into different cohorts for the same period. The trick is to not arrange your age / period / cohorts so that the linear relationship is preserved. Commented Jan 19, 2023 at 2:50

I might just add the features and fit on that one feature instead of both. That forces the coefficient on each to be the same, and you can run your usual software functions for a Poisson regression.

$$y=\beta_0+\beta_1x_1 +\beta_2x_2\\ \Updownarrow\\ y=\beta_0+\beta_1\left( x_1+x_2 \right)$$

• Thanks you! It is ok to do that? I want to use the glm to do prediction.Although I had thought of this method, I am very confused whether I should use or not. Commented Jan 18, 2023 at 3:10
• Also, may I ask one more thing: Will the confidence intervals of this two variables be the same...? Commented Jan 18, 2023 at 3:13
• @gczday There’s only one parameter, so I am struggling to understand the question. Do you want different confidence intervals for the common parameter on each variable? Why?
– Dave
Commented Jan 18, 2023 at 3:14
• Since I am going to plot a graph of relative risk of period effects, where confidence intervals are needed for reference group and the group that was assumed to have same linear effect as reference group Commented Jan 18, 2023 at 3:17
• It doesn’t make sense to me why there would be two confidence intervals, but maybe I’m stuck thinking of the constrained regression with two parameters as being the same as the regression on a feature that is the sum of the original features. (Actually, those sound exactly the same, though I confess that I have never considered this before now.) // You might want to edit your original question to include your interest in the confidence interval and if it makes sense to have two.
– Dave
Commented Jan 18, 2023 at 3:30