4
$\begingroup$

I am fitting a poisson GLM to model claim rates. Since I have 1.5m+ records, I have aggregated my data (to improve efficiency). My understanding is that using aggregated data with a poisson GLM will not effect the coefficient estimates. Indeed, if I fit a GLM with no regularization (and using the log of exposure as an offset) I can exactly recreate the same coefficient estimates using either the full dataset or aggregated data.

However, as soon as I introduce a lambda penalty, the coefficients change depending on whether I use full individual or aggregated data. The larger the penalty, the larger the difference in coefficients.

My question is - should I be able to get the same coefficient estimates using aggregated data with a lambda penalty as I do with individual data? And finally, if not, is it 'ok' to use regularization with aggregated data?

For info, i'm using the GLM function from H2O.ai, in R Studio. Here is some reproducible code (for simplicity i've only simulated 50k records here):

library(simstudy)
set.seed(123)

# simulate a dataset with 2 variables - 'Male' and 'Target'
def = defData(varname = "Male", dist = "binary", formula = 0.4)
def = defData(def, varname = "Target",dist = "binary",formula = "ifelse(Male == 1, 0.3, 0.7)")

# generate the dataset
dt = genData(50000,def)

# variable to indicate that each row represents 1 exposure period
dt$n = 1

dt$Male = as.factor(dt$Male)
dt = as.data.frame(dt)

# create an aggregated version of data
dt2 = aggregate(x = dt[c('Target','n')],
                by = list(Male = dt$Male),
                FUN = sum)

# add an offset column to aggregated data - log of exposure
dt2$offset = log(dt2$n)

# initialise h2o
library(h2o)
h2o.init(nthreads = -1)
dt = as.h2o(dt)
dt2 = as.h2o(dt2)

# fit GLM to full data, no regularization
mod1 = h2o.glm(x = 'Male', y = 'Target',training_frame = dt,family = 'poisson',
               lambda = 0,alpha = 0,seed = 123)
round(mod1@model$coefficients,5)

#Intercept    Male.1 
#-0.35798  -0.83271 

# mod2 - use aggregated data with an offset
mod2 = h2o.glm(x = 'Male', y = 'Target',training_frame = dt2,family = 'poisson',
               lambda = 0,alpha = 0,offset_column = 'offset',seed = 123)
round(mod2@model$coefficients,5)

#Intercept    Male.1 
#-0.35798  -0.83271 

# now repeat, with regularization. First, full data
mod3 = h2o.glm(x = 'Male', y = 'Target',training_frame = dt,family = 'poisson',
               lambda = 0.01,alpha = 0,seed = 123)
round(mod3@model$coefficients,5)

#Intercept    Male.0    Male.1 
#-0.76286   0.39547  -0.39547 

# now on aggregated data
mod4 = h2o.glm(x = 'Male', y = 'Target',training_frame = dt2,family = 'poisson',
               lambda = 0.01,alpha = 0,offset_column = 'offset',seed = 123)
round(mod4@model$coefficients,5)

#Intercept    Male.0    Male.1 
#-0.77433   0.41635  -0.41635

UPDATE:

I have been through the H2O GLM documentation and can't find anything obviously helpful in this case. I have also tried removing collinear columns in the models above, removing any standardizing and ensuring that the number of training iterations is the same for each version of the model - all with no luck in terms of reproducibility.

Here is a link to the h2o GLM documentation: http://docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/glm.html

In addition, I have recreated the same issue using the R 'glmnet' package as well as h2o. I can post code if necessary.

(Have amended question to reflect that i'm modelling rate data)

UPDATE:

I'm also recreated the issue using the glmnet R package. See below

 # using a different GLM package - glmnet
    library(simstudy)
    # simulate a dataset
    set.seed(123)

# simulate a dataset with 3 variables - 'Male', 'Occupation' and 'Target'
def = defData(varname = "Male", dist = "binary", formula = 0.4)
def = defData(def,varname = "Occupation", dist = "categorical", formula = "0.2;0.4;0.4")
def = defData(def, varname = "Target",dist = "binary",formula = "ifelse(Male == 1, 0.3, 0.7)")

# generate the dataset
dt = genData(50000,def)

# variable to indicate exposure for each record
n = runif(nrow(dt), 1,3)
dt$n = n

dt$Male = as.factor(dt$Male)
dt$Occupation = as.factor(dt$Occupation)
dt = as.data.frame(dt)

# create an aggregated version of data
dt2 = aggregate(x = dt[c('Target','n')],
                by = list(Male = dt$Male,
                          Occupation = dt$Occupation),
                FUN = sum)

library(glmnet)

x1 = as.matrix(dt[2:3])
y1 = as.matrix(dt$Target)

fit1 = glmnet(x = x1,y = y1,family = "poisson",alpha = 0, lambda = 0.1, offset = log(dt$n), standardize = F)
coef(fit1)

# now using aggregated data
x2 = as.matrix(dt2[1:2])
y2 = as.matrix(dt2$Target)

fit2 = glmnet(x = x2,y = y2,family = "poisson",alpha = 0, lambda = 0.1, offset = log(dt2$n), standardize = F)
coef(fit2)

UPDATE:

I have been doing a bit more digging into this for anyone who's interested. I used the 'lambda_search' option set to TRUE (in h2o.glm) to see how h2o would automatically search over lambda to find the best model fit. I then used h2o.getGLMFullRegularizationPath() to check the different values of lambda that h2o was testing. I noticed that the lambdas h2o selected for the aggregated model were always 25000 times the value of the lambdas selected for the original model. Since my dataset has been reduced by a factor of 25000 (50000 records to 2 records) it seems my lambda needs to change accordingly. Indeed, when i've tested the two models again I can recreate the regularized original model almost exactly using aggregated data if lambda is set to be original lambda * nrow(full data)/nrow(aggregated data). For example:

mod3 = h2o.glm(x = 'Male', y = 'Target',training_frame = dt,family = 'poisson',
               lambda = 0.1,alpha = 0,seed = 123, offset_column = 'offset',standardize = F)
mod3@model$coefficients,5

# Gives result:
#Intercept     Male.0     Male.1 
#-1.3984380  0.2773347 -0.2773347

Then, if I try and recreate this model using aggregated data with regularization:

# now on aggregated data
lambda = 0.1 * nrow(dt)/nrow(dt2)
mod4 = h2o.glm(x = 'Male', y = 'Target',training_frame = dt2,family = 'poisson',
               lambda = lambda,alpha = 0,offset_column = 'offset',seed = 123,standardize = F)
mod4@model$coefficients

# gives result:
# Intercept     Male.0     Male.1 
# -1.3984384  0.2773349 -0.2773349

I am not quite sure how the maths works here so if anyone could help explain that it would be much appreciated.

Improve this question
$\endgroup$
5
  • $\begingroup$ There should be no problem in principle, but in practice it depends on the implementation ... $\endgroup$
    – kjetil b halvorsen ♦
    Feb 19, 2018 at 11:02
  • $\begingroup$ @kjetilbhalvorsen thanks - so do you mean it's likely to be something related to the R package i'm using to fit the GLM? docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/glm.html $\endgroup$
    – Joshua Olney
    Feb 19, 2018 at 12:05
  • $\begingroup$ Yes, probably. Which package are you using? $\endgroup$
    – kjetil b halvorsen ♦
    Feb 19, 2018 at 12:29
  • $\begingroup$ @kjetilbhalvorsen I'm using the GLM function in the 'h2o' R package. I could post reproducible code if it helps - or move to Stack Overflow if it might be a programming issue? thanks $\endgroup$
    – Joshua Olney
    Feb 19, 2018 at 13:42
  • 1
    $\begingroup$ Your last update is interesting, good work, and you have very probably found the culprit. That is sometimes done to make the two terms of the objective function comparable (the first term, sum of squares (or deviance) will obviously increase with $N$, while the second term $\lambda$ times the norm of $\beta$, does not depend on $N$, so to make this comparable over different datasets with different sample size, the last term is multiplied by $N$. I will update my answer. $\endgroup$
    – kjetil b halvorsen ♦
    Feb 21, 2018 at 17:30

1 Answer 1

Reset to default
2
$\begingroup$

I looked at this, and have some comments. First, you use offsets to represent weights. This are not necessarily the same, so I changed to use weights_column argument. See Difference between: Offset and Weights? (my modified version of your code is given below). But there are other problems. Two, I included the argument earl_stopping=FALSE just to be sure of comparability. Third, there is a standardize argument with default value TRUE, it should be set to FALSE. If you want to standardize it must be done BEFORE aggregation, else the coefficients will have different meaning in the two cases. But still, I get different results, to me it seems like the standardize argument have no effect. That would be an error in the h20 package which you should investigate and report. And maybe try some other package, like glmnet.

My modified code:

library(simstudy)  # On CRAN 
set.seed(7*11*13)

# simulate a dataset with 2 variables - 'Male' and 'Target'
def <- defData(varname = "Male", dist = "binary", formula = 0.4)
def <- defData(def, varname = "Target",dist = "binary",formula = "ifelse(Male == 1, 0.3, 0.7)")

# generate the dataset
dt <- genData(50000,def)

# variable to indicate that each row is 1 life
dt$n = 1

dt$Male <- as.factor(dt$Male)
dt <- as.data.frame(dt)

# create an aggregated version of data
dt2 <- aggregate(x = dt[c('Target','n')],
                by = list(Male = dt$Male),
                FUN = sum)

# add a weight column to aggregated data  exposure
dt2$weight <- dt2$n

# initialise h2o
library(h2o)  on CRAN  
h2o.init(nthreads = -1)
dt <-  as.h2o(dt)
dt2 <- as.h2o(dt2)

# fit GLM to full data, no regularization
mod1 <- h2o.glm(x = 'Male', y = 'Target',training_frame = dt,family = 'poisson', lambda = 0,alpha = 0,seed = 7*11*13, early_stopping=FALSE, standardize=FALSE)
round(mod1@model$coefficients,5)

# mod2 - use aggregated data with a weight column
mod2 <- h2o.glm(x = 'Male', y = 'Target', training_frame = dt2, family = 'poisson', lambda = 0, alpha = 0, weights_column = 'weight', seed = 7*11*13, early_stopping=FALSE, standardize=TRUE)
round(mod2@model$coefficients,5)

# now repeat, with regularization. First, full data
mod3 <- h2o.glm(x = 'Male', y = 'Target',training_frame = dt,family = 'poisson',
               lambda = 0.01,alpha = 0,seed = 7*9*13, early_stopping=FALSE, standardize=FALSE)
round(mod3@model$coefficients,5)

#Intercept    Male.0    Male.1 
#-0.76286   0.39547  -0.39547 

# now on aggregated data
mod4 <- h2o.glm(x = 'Male', y = 'Target',training_frame = dt2,family = 'poisson', lambda = 0.01,alpha = 0, weights_column = 'weight', seed = 7*11*13, early_stopping=FALSE, standardize=FALSE)
round(mod4@model$coefficients,5)
Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ h2o package has h2o.coef and h2o.coef_norm to get the coefficient where former returns the coefficient before standardisation. $\endgroup$
    – Sixiang.Hu
    Feb 20, 2018 at 13:43
  • $\begingroup$ perhaps I should have been clearer as to why I was using an offset - where I've said 'add a variable to indicate that each row is 1 life' what i really meant is that in the full dataset, each row represents someone with exposure of period 1 (e.g 1 year) - but equally this could have been any number. I then didn't include this explicitly as an offset in mod1. My understanding is then that I want to include the log of exposure when using my aggregated data as I want the model to know that my 'Target' varies as exposure time varies? I suppose then that i'm really modelling a rate? $\endgroup$
    – Joshua Olney
    Feb 20, 2018 at 14:09
  • $\begingroup$ Yes, then you are modelling the rate, so what you did looks correct then. But still what I said about standardization is correct, so you should investigate that. $\endgroup$
    – kjetil b halvorsen ♦
    Feb 20, 2018 at 14:16
  • 1
    $\begingroup$ Yes thanks - I'll flag that on the h2o google group as the standardize condition seems to make no difference to me also $\endgroup$
    – Joshua Olney
    Feb 20, 2018 at 14:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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