[R] Negative Binomial Regression
Ben Bolker
bolker at ufl.edu
Tue Jul 29 18:11:26 CEST 2008
BEGIN PGP SIGNED MESSAGE
Hash: SHA1
Prof Brian Ripley wrote:
 On Tue, 29 Jul 2008, Ben Bolker wrote:

> jcarmichael <jcarmichael314 <at> gmail.com> writes:
>
>>
>>
>> Hello.
>>
>> I am attempting to duplicate a negative binomial regression in R.
>> SAS uses
>> generalized estimating equations for model fitting in the GENMOD
>> procedure.
>>
>> proc genmod data=mydata (where=(gender='F'));
>> by agegroup;
>> class id gender type;
>> model count = var1 var2 var3 /dist=NB link=log offset=lregtm;
>> repeated subject=id /type=exch;
>> run;
>>
>> Since my dataset has several observations for each subject, I need the
>> REPEATED statement in order to indicate dependence among observations
>> with
>> the same subject ID and independence amongst those with distinct subject
>> IDs. The TYPE statement goes on to specify the structure of the
>> correlation
>> matrix to be used (exchangeable in this case).
>
> I would try glmmPQL in the MASS package. I don't think you
> can *quite* get negative binomial regression this way, but
> you can definitely get a quasipoisson model. I think exchangeable
> correlation corresponds to correlation=corCompSymm() in your
> glmmPQL command.

 The problem here is that GLMM and GEE are not fitting the same model 
 in one the coefficients are subjectspecific and in the other
 populationaverage (see MASS4 or Diggle, Liang, Zeger +/ Heagarty).

 There are several R packages for GEE, including gee, yags, geepack. The
 documentation of geeglm (geepack) claims it can be used with families as
 in glm(), so you could try it with MASS's negative.binomial family.

~ Point taken (although I guess I was pointing the original poster
to a way to do a reasonable analysis, not necessarily to duplicate
the SAS analysis as requested). Will the negative.binomial family
really work for this, since it seems to require a fixed theta
(overdispersion) parameter?
~ If I very naively do the following:
library(geepack)
data(dietox)
mf2 < formula(Weight~Cu*Time+I(Time^2)+I(Time^3))
gee2 < geeglm(mf2, data=dietox, id=Pig,
~ family=poisson("identity"),corstr="ar1")
library(MASS)
gee2 <
geeglm(mf2,data=dietox,id=Pig,family=negative.binomial(theta=100),corstr="ar1")
~ gives an error "variance invalid" 
~ so the whole thing would seem to take a bit of troubleshooting
~ (geeglm also gives warnings about noninteger Poisson values  I
don't know why a Poisson link is being used in this example for
a noninteger Weight value ... ?)
~ Ben Bolker
BEGIN PGP SIGNATURE
Version: GnuPG v1.4.6 (GNU/Linux)
Comment: Using GnuPG with Mozilla  http://enigmail.mozdev.org
iD8DBQFIj0Euc5UpGjwzenMRAhs/AJ9EOu8RpJibN9LMLbI0p1tCTls5xQCfRLnf
6hkinpQu0PMQgh4f/Sl5BaM=
=y7CX
END PGP SIGNATURE
More information about the Rhelp
mailing list