R not Matching Stata

irkgreen

New Member
Hello,

I am trying to perform a negative binomial regression in RStudio on a dataset but I want to specify a formula for overdispersion (shape parameter). I am able to match almost exactly what Stata produces but I am getting the negative of what Stata produces.

My negative binomial model form is: y=a*ADT^b+L^c and the shape parameter is theta=const*L^b1

I don't think I can use glm.nb because it doesn't have an option for a variable dispersion (as far as I know).

My R code is:
Code:
library(MASS)
library(gnlm)

CSVpath = "MYPATH/Urb2Undiv_D2.csv"

#Point to variables
y <- c(data$Total) lnADT <- c(log(data$ADT))
lnL <- c(log(data$Length)) Lx <- c(data$Length)

SPF = gnlr(y, dist="negative binomial", mu=~exp(a+b*lnADT+c*lnL), shape=~(const+b1*lnL), pmu=list(a=-1,b=.4,c=.7), pshape=c(0,0))
Notice the coefficients for lnAlpha (dispersion): in R const=1.2440 and in Stata cons=-1.244108

Is this a bug in R (Stata's results seem more intuitive)?

My R output is:

Call:
gnlr(y, dist = "negative binomial", mu = ~exp(a + b * lnADT +
c * lnL), shape = ~(const + b1 * lnL), pmu = list(a = -1,
b = 0.4, c = 0.7), pshape = c(0, 0))

negative binomial distribution

Response: y

Log likelihood function:
{
m <- mu1(p)
t <- sh1(p)
s <- exp(t)
-sum(wt * (lgamma(y + s) - lgamma(s) + s * t + y * log(m) -
(y + s) * log(s + m)))
}

Location function:
~exp(a + b * lnADT + c * lnL)

Log shape function:
~(const + b1 * lnL)

-Log likelihood 3910.424
Degrees of freedom 1933
AIC 3915.424
Iterations 34

Location parameters:
estimate se
a -1.0433 0.21282
b 0.3638 0.02558
c 0.6717 0.04123

Shape parameters:
estimate se
const 1.2440 0.2091
b1 0.3937 0.1216

Correlations:
1 2 3 4 5
1 1.00000 -0.94912 0.09922 0.04912 0.05547
2 -0.94912 1.00000 0.20187 -0.03999 -0.04530
3 0.09922 0.20187 1.00000 0.03079 0.03404
4 0.04912 -0.03999 0.03079 1.00000 0.95371
5 0.05547 -0.04530 0.03404 0.95371 1.00000

>

My Stata Output is:

Generalized negative binomial regression Number of obs = 1938
Wald chi2(2) = 400.28
Prob > chi2 = 0.0000
Log pseudolikelihood = -3910.4238 Pseudo R2 = 0.0432