[ADMB Users] Poisson GLM with a linear link function

[Mollie Brooks]

Hi Mark,
I was able to do a model like you're talking about (Poisson with identity link) without constraining the mean. I also had random effects of year and individual. I had to start it with good values. You could first try a linear model (with a normal distribution). Check if the linear model ever goes below zero. In my case, the mean was always positive.

For profiling I did have to constrain the mean using posfun (see code at end of message). Some purists would say that if you have to use posfun, then you aren't actually fitting the statistical model you say you are. I'm pragmatic about it, but I now people should be cautious. You should make the program notify you if posfun is having to make alterations by checking if fpen>0. Then think seriously about the alterations.

After doing that, I learned that to be kosher I could have followed the method in this paper. Check it out:
Marschner, I. C. (2010). Stable Computation of Maximum Likelihood Estimates in Identity Link Poisson Regression. Journal of Computational and Graphical Statistics, 19(3), 666–683. doi:10.1198/jcgs.2010.09127

At least I can point to this paper (and others) to say that Poisson models with identity links are legit.

cheers,
Mollie

DATA_SECTION
	init_int nobs;
	init_int nyears;
	init_vector shts_next(1,nobs);
	init_vector shts(1,nobs);
	init_int nfix;
	init_matrix fixmm(1,nobs,1,nfix);
	init_int nind;
	init_ivector indiv(1,nobs);
	init_ivector year(1,nobs);
	
PARAMETER_SECTION
	init_vector fcoeffs(1,nfix,1);
	init_bounded_number log_indiv_std(-25,10,1);
	init_bounded_number log_year_std(-25,10,3);
	sdreport_number indiv_std1;
	random_effects_vector rcoeffs_indiv(1,nind,1);
	random_effects_vector rcoeffs_year(1,nyears,3);
	objective_function_value jnll;

PROCEDURE_SECTION
	jnll=0.0;
	for(int i=1; i<=nobs; i++)
	{
		obs(i,fcoeffs,log_indiv_std,log_year_std,rcoeffs_indiv(indiv(i)), rcoeffs_year(year(i)));
	}
	//now add in the nll of the random effects
	for(int i=1; i<=nind; i++)
	{
		rand(rcoeffs_indiv(i));
	}
	for(int i=1; i<=nyears; i++)
	{
		rand(rcoeffs_year(i));
	}
	indiv_std1=exp(log_indiv_std);
	
SEPARABLE_FUNCTION void obs(int i,const dvar_vector& fcoeffs,const dvariable& log_indiv_std,const dvariable& log_year_std,const dvariable& rcoeff_indiv,const dvariable& rcoeff_year)
	dvariable indiv_std = exp(log_indiv_std);
	dvariable year_std = exp(log_year_std);
	dvariable prediction=fixmm(i)*fcoeffs+
	 					indiv_std*rcoeff_indiv+
						year_std*rcoeff_year;
	dvariable fpen=0.0;
	
	prediction=1.e-8+posfun(prediction,.01,fpen);
	
	jnll+=fpen;

	jnll+=-shts_next(i)*log(prediction)+prediction + gammln(shts_next(i)+1.);
	
SEPARABLE_FUNCTION void rand(const dvariable& ui)
	jnll+=0.5*(ui*ui)+0.5*log(2.*M_PI);
	
TOP_OF_MAIN_SECTION
	arrmblsize = 30000000;
    gradient_structure::set_MAX_NVAR_OFFSET(168514);


Mollie Brooks
Ph.D. Candidate
NSF IGERT Fellow
Biology Department
University of Florida
mbrooks at ufl.edu
http://people.biology.ufl.edu/mbrooks




On 19 Jun 2012, at 11:41 AM, Mark Payne wrote:

> Dear ADMBers,
> 
> I have a problem that is essentially a poisson GLM with a linear link function i.e.
> 
> y | lambda ~ Pois(lambda)
> 
> lamba = a*x1 + b*x2 + c*x3.....
> 
> where  y are my observations, lambda is the expected value, x are the explanatory variables, and a-c are the coefficients that I want to estimate. Now, a linear link function is relatively unusual in this context, but it arises from the process that I am modelling and therefore is a fixed part of the model - essentially the underlying process giving rise to the observations is a linear sum by nature, and I am then observing the results via a counting process.
> 
> Unfortunately, the problem is that the mean value in the poisson likelihood, lambda, needs to be kept strictly positive, otherwise the likelihood is undefined. The question is, how can I best do this in ADMB *and* maintain the linear structure of my model? Are there any tricks in ADMB or transforms that could be used to make this problem more stable?
> 
> Best wishes,
> 
> Mark
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20120619/c1ae6e0f/attachment.html>