r-spatial / spatialreg

spatialreg: spatial models estimation and testing
https://r-spatial.github.io/spatialreg/
43 stars 11 forks source link

Using the ME function for count data of negative binominal distribution #34

Closed zhenchun closed 1 year ago

zhenchun commented 1 year ago

I have some data here available at: https://github.com/zhenchun/data

require(spdep)
require(spatialreg)
require(sf)

nc = st_read("nc.shp", package="sf")

nb <- poly2nb(nc, queen=TRUE)

wt<-nb2listw(nb,zero.policy=TRUE, style="W")

fit<-glm.nb(n~pnhw13_17+offset(log(tt13_17)), nc, na.action=na.omit)
> summary(fit)

#Call:
#glm.nb(formula = n ~ pnhw13_17 + offset(log(tt13_17)), data = nc, 
#    na.action = na.omit, init.theta = 0.6835949402, link = log)

#Deviance Residuals: 
#    Min       1Q   Median       3Q      Max  
#-2.4175  -1.1024  -0.5685   0.1666   4.1108  

#Coefficients:
 #           Estimate Std. Error z value Pr(>|z|)    
#(Intercept)  -8.2594     0.1569  -52.63   <2e-16 ***
#pnhw13_17     3.0596     0.2121   14.43   <2e-16 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#(Dispersion parameter for Negative Binomial(0.6836) family taken to be 1)

#    Null deviance: 1132.94  on 800  degrees of freedom
#Residual deviance:  910.99  on 799  degrees of freedom
#AIC: 5665.7

#Number of Fisher Scoring iterations: 1

#              Theta:  0.6836 
#          Std. Err.:  0.0390 

# 2 x log-likelihood:  -5659.6900 

me.fit<-spatialreg::ME(n~pnhw13_17,offset=offset(log(tt13_17)), data = nc, family=negative.binomial(0.68236, link=log),zero.policy=TRUE, listw = wt, verbose=T,alpha=.05,  na.action=na.omit)

#Error in glm.fit(x = iX, y = Y, weights = weights, offset = offset, family = family) : 
 # NA/NaN/Inf in 'x'

Anybody knows how this error come from? Thanks a lot for your help.

I also learnt from these links https://www.mail-archive.com/r-sig-geo@r-project.org/msg03967.html https://rpubs.com/corey_sparks/111362

rsbivand commented 1 year ago

Yes, ME() does not support the Negative Family. The mailing list posting you found includes a modified version of ME(), but nothing was ever tested or added to the original function. Please see spmoran https://cran.r-project.org/package=spmoran as an actively maintained package for Moran eigenvectors.

zhenchun commented 1 year ago

Thanks a lot for your help!