[R] Constrained Poisson model / Bayesian Poisson model

[David Winsemius]

> On Jan 23, 2016, at 5:24 AM, mara.pfleiderer at uni-ulm.de wrote:
> 
> Hi David,
> 
> I'm sorry. I'm not familiar with posting problems on helppages.
> As the data I deal with is confidential I can't provide all details,
> but I will try to be as precise as possible about my problem:
> 
> As I said I'm working on a Poisson regression model with a linear
> predictor and the identity as link function.

Then that is not what most people would call a "Poisson regression model".

In particular: I have data which consists of observations of two random variables X and Y over the period of about 3 months. I have clustered this data into hours, so I received 2088 observations of those variables and each observation represents the values of the variables in one hour. Now I assume Y to be Poisson distributed (Y_i~Poi(lambda_i)) and that the values of Y come from two different impacts, so I split Y_i into two Poisson-distributed variables Y_i1 and Y_i2. I assume that the values of X have a long-term effect (of at least one day) on the part Y_i1 and I have estimators for the parameters lambda_i2, but I have no exact values of the variables of Y_i1 and Y_i2 but only of Y_i as a whole.

> So my model looks as follows:
> fml=Y_i1 ~ b_1*X_i+....+b_n*X_(i-n+1) - lambda_i2 -1 with n>=24
> That means that the last 24 (or more) values of X influence the value of Y_i1 in an additive way and I have no intercept(b_0=0).
> Now I made a matrix whose rows represent Y_i, all of its 24 (or more) regressors for each observation of Y_i and the corresponding estimator lambda_i2.
> Then I used glm(fml, family=poisson(link="identity"), data=matrix) and tried it for different values of n (=24,48,36,...).

I worry that might not construct the model you described above. The use of `poisson` constructs Poisson distributed _errors_ with an additive link, and to my understanding does _not_ assume Poisson distributed Y values.


> But always some of the coefficients received negative values which doesnt't make sense in interpretation. (The values of X represent certain events which can only have a positive or none effect on the value of Y.)
> 
> Now I want to use the constraint b_i >=0 in my model, but I don't know how I can do this.
> 
> I also thought of analyzing this model in a Bayesian way, but yet I haven't found a Bayesian version of glm() for the Poisson distribution such that I can specify the prior for the b_i on my own. (Then I could include the positivity in the prior.) Do you have a hint for me?

The R blogoshere has been heralding the arrival of the `rstanarm` package which provides a glm-like interface to an MCMC (Bayesian) estimation engine. Version 2.90 is on CRAN. The list of authors in the description file is impressive: 

Authors at R: c(person("Jonah", "Gabry", email = "jsg2201 at columbia.edu", role = "aut"),
             person("Trustees of", "Columbia University", role = "cph"),
             person("R Core", "Deveopment Team", role = "cph", 
                     comment = "R/pp_data.R, R/stan_aov.R"),
             person("Douglas", "Bates", role = "cph", comment = "R/pp_data.R"),
             person("Martin", "Maechler", role = "cph", comment = "R/pp_data.R"),
             person("Ben", "Bolker", role = "cph", comment = "R/pp_data.R"),
             person("Steve", "Walker", role = "cph", comment = "R/pp_data.R"),
             person("Brian", "Ripley", role = "cph", 
                    comment = "R/stan_aov.R, R/stan_polr.R"),
             person("William", "Venables", role = "cph", comment = "R/stan_polr.R"),
             person("Ben", "Goodrich", email = "benjamin.goodrich at columbia.edu", 
                    role = c("cre", "aut")))

> 
> I'm sorry if I went too much into detail now...
> I hope you understand my point now and have some answers for me!
> 
> Best,
> Mara
> 
> 
> Zitat von David Winsemius <dwinsemius at comcast.net>:
> 
>> 
>>> On Jan 22, 2016, at 7:01 AM, mara.pfleiderer at uni-ulm.de wrote:
>>> 
>>> Hi all,
>>> 
>>> I am dealing with a problem about my linear Poisson regression   model (link function=identity).
>>> 
>>> I am using the glm()-function which results in negative   coefficients, but a negative influence of the regressors wouldn't   make sense.
>> 
>> Negative coefficients merely indicate a lower relative rate. You   need to be more specific about the exactly data and model output   before you can raise our concern to a level where further comment   can be made.
>> 
>> 
>>> 
>>> (i) Is there a possibility to set constraints on the regression   parameters in glm() such that all coefficients are positive? Or is   there another function in R for which this is possible?
>>> 
>>> (ii) Is there a Bayesian version of the glm()-function where I can   specify the prior distribution for my regression parameters? (e.g.   a Dirichlet prior s.t. the parameters are positive)
>>> 
>>> All this with respect to the linear Poisson model...
>> 
>> As I implied above, the word "linear" means something different than   "additive" when the link is log().
>> 
>> --
>> 
>> David Winsemius
>> Alameda, CA, USA
>> 

David Winsemius
Alameda, CA, USA