3
$\begingroup$

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:

enter image description here

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?

$\endgroup$
5
  • 1
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Commented Jan 18, 2023 at 16:42
  • $\begingroup$ Thank you for pointing out this! I am not familiar with how to construct sample data... $\endgroup$
    – doraemon
    Commented Jan 19, 2023 at 2:12
  • 1
    $\begingroup$ You will see an example in my updated answer below. $\endgroup$
    – jbowman
    Commented Jan 19, 2023 at 2:34
  • 2
    $\begingroup$ 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!) $\endgroup$ Commented Feb 27 at 3:38
  • $\begingroup$ Thank you for your reminder! I have asked a new question. $\endgroup$
    – doraemon
    Commented Feb 27 at 3:49

2 Answers 2

7
$\begingroup$

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.

Expanded answer in response to comments:

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.

$\endgroup$
14
  • $\begingroup$ 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. >~< $\endgroup$
    – doraemon
    Commented Jan 18, 2023 at 5:01
  • 2
    $\begingroup$ 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." $\endgroup$
    – jbowman
    Commented Jan 19, 2023 at 1:56
  • 1
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Commented Jan 19, 2023 at 2:06
  • 1
    $\begingroup$ I've expanded my answer with an example above! $\endgroup$
    – jbowman
    Commented Jan 19, 2023 at 2:28
  • 1
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Commented Jan 19, 2023 at 2:50
6
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – doraemon
    Commented Jan 18, 2023 at 3:10
  • $\begingroup$ Also, may I ask one more thing: Will the confidence intervals of this two variables be the same...? $\endgroup$
    – doraemon
    Commented Jan 18, 2023 at 3:13
  • 1
    $\begingroup$ @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? $\endgroup$
    – Dave
    Commented Jan 18, 2023 at 3:14
  • $\begingroup$ 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 $\endgroup$
    – doraemon
    Commented Jan 18, 2023 at 3:17
  • 2
    $\begingroup$ 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. $\endgroup$
    – Dave
    Commented Jan 18, 2023 at 3:30

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.