2
$\begingroup$

I'm currently using bayesian logistic regression for binary clasification (arm package - R) , one of its many wonderful features that I love to exploit is the possibility the update my coefficient every time a new batch of data arrives (What I understand as bayesian updating), avoiding the need to upload all my data at once.

I'm running a sanity check to make sure my code works properly. To my understanding, running two blocks of bayesian logistic regressions, where the first block's coefficients will be used as priors to the second block, will result in a similar output as if I would have run all my data at once (or am I wrong !?)

Sanity check code:

library(data.table)
library(arm)

#some sample data:
n <- 100
x1 <- rnorm (n)
x2 <- rbinom (n, 1, .5)
b0 <- 1; b1 <- 1.5; b2 <- 2
y <- rbinom (n, 1, invlogit(b0+b1*x1+b2*x2))
z1 <- trunc(runif(n, 4, 9))
dt<-data.table(y=y,x1=x1,z1=z1)
head(dt)

#i like to work in agregated format :)
first_batch<-dt[,.(count=.N),.(y,z1)]

#full data - without beysian updating 
full_dt<-rbind(first_batch,first_batch)
full_dt<-full_dt[,.(count=sum(count)),.(y,z1)]
full_model<- bayesglm(y ~ factor(z1),weights=count,family=binomial(link="logit"),data=full_dt)
#The following is what i really want to compute:

#using first batch
prior_estimates <-bayesglm(y ~ factor(z1),weights=count,family=binomial(link="logit"),data=first_batch)
prior_estimates_dt<-data.table(term=attr(summary(prior_estimates)$coefficients,"dimnames")[[1]],
                         estimate=as.data.table(summary(prior_estimates)$coefficients)$Estimate,
                         std.error=as.data.table(summary(prior_estimates)$coefficients)[,2,with = FALSE])
postirior_model <-bayesglm(y ~ factor(z1),
                   prior.mean=prior_estimates_dt[term!="(Intercept)"]$estimate,
                   prior.scale=prior_estimates_dt[term!="(Intercept)"]$std.error,
                   weights=count,family=binomial(link="logit"),data=first_batch)

#compate resutls
display (full_model,detail=FALSE)
display(prior_estimates,detail=FALSE)
display(postirior_model,detail=FALSE)

> display (full_model,detail=FALSE)
bayesglm(formula = y ~ factor(z1), family = binomial(link = "logit"), 
    data = full_dt, weights = count)
            coef.est coef.se
(Intercept)  0.72     0.31  
factor(z1)5  1.67     0.71  
factor(z1)6 -0.30     0.42  
factor(z1)7  0.30     0.47  
factor(z1)8 -0.42     0.43  
---
n = 10, k = 5
residual deviance = 236.4, null deviance = 250.7 (difference = 14.3)
> display(prior_estimates,detail=FALSE)
bayesglm(formula = y ~ factor(z1), family = binomial(link = "logit"), 
    data = first_batch, weights = count)
            coef.est coef.se
(Intercept)  0.73     0.42  
factor(z1)5  1.52     0.91  
factor(z1)6 -0.31     0.57  
factor(z1)7  0.27     0.63  
factor(z1)8 -0.42     0.59  
---
n = 10, k = 5
residual deviance = 118.3, null deviance = 125.4 (difference = 7.1)
> display(postirior_model,detail=FALSE)
bayesglm(formula = y ~ factor(z1), family = binomial(link = "logit"), 
    data = first_batch, weights = count, prior.mean = prior_estimates_dt[term != 
        "(Intercept)"]$estimate, prior.scale = prior_estimates_dt[term != 
        "(Intercept)"]$std.error)
            coef.est coef.se
(Intercept)  0.73     0.29  
factor(z1)5  1.64     0.62  
factor(z1)6 -0.32     0.35  
factor(z1)7  0.29     0.40  
factor(z1)8 -0.43     0.36  
---
n = 10, k = 5
residual deviance = 118.2, null deviance = 125.4 (difference = 7.2)

Notice the results of full_model VS. postirior_model , the bayesian update seemed to get the posterior results closer to the full models results, but not exactly as I anticipated, can anyone explain the source of this mismatch?

Thanks in advance!

$\endgroup$

1 Answer 1

1
$\begingroup$

My guess is that when you are performing the Bayesian updating, you are not using the correct updated prior. Logistic regression has no conjugate prior thus you will not be able to perform Bayesian updating the way that you want to. In addition, the default prior for bayesglm is a Cauchy distribution and you have only updated the location and scale for this distribution (in fact you haven't ever updated the intercept location and scale), but have never updated the degrees of freedom. Finally, the prior in bayesglm is always independent amongst the coefficients, but the updated posterior will not necessarily be independent. (There is also the question of what the scaled argument is doing to the prior, so you will likely want scaled=FALSE.)

You may be able to get closer to what you want by taking samples from the posterior and then optimizing the location, scale, and degrees of freedom to find the best fitting t distribution (likely you will do this for each coefficient separately). Then use this as the prior for the next data subset.

$\endgroup$
1
  • $\begingroup$ thanks for this great answer, just a short clarification, assuming ill take care if the degrees of freedom, will my only concern is that my posterior coefficient are independent? because in my real data is seems to be the case. $\endgroup$ Sep 1, 2016 at 8:54

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.