est.sc: Fitting spatial count regression models
In spatcounts: Spatial count regression

Description Usage Arguments Value References See Also Examples

MCMC algorithm for the Poisson, the GP, the NB, the ZIP and the ZIGP regression models with spatial random effects.

est.sc(Yin, fm.X, region, model = "Poi", gmat, nmat, totalit, 
fm.ga = TRUE, t.i = NULL, phi0 = 1, omega0 = 0, r0 = 1, 
beta0 = NULL, gamma0 = NULL, sigma0 = 1, psi0 = 1, Tau = 10, 
alpha = 2)

`Yin`	response vector of length n.
`fm.X`	formula for mean design.
`region`	region of each observation as vector of length n.
`model`	the regression model. Currently, the regression models "Poi", "NB", "GP", "ZIP" and "ZIGP" are supported. Defaults to 'Poi'.
`gmat`	spatial adjacency matrix, where entry (i,j) is 1 if region i is a neighbor of region j and 0 else. See data(sim.gmat) for an example.
`nmat`	matrix containing the number of neighbors of each region (last column) and the neighbors of each region (first (maximual number of neighbours) columns), filled up by zeros. See data(sim.nmat) for an example.
`totalit`	number of MCMC iterations, i.e. length of the Markov chain.
`fm.ga`	should the spatial random effects be included (defaults to TRUE)?
`t.i`	exposure vector.
`phi0`	starting value for the over-dispersion parameter for GP and ZIGP model.
`omega0`	starting value for the extra proportion for ZIP and ZIGP model.
`r0`	starting value for the scale paramter for NB model.
`beta0`	starting values for the regression parameters.
`gamma0`	starting values for the spatial paramters.
`sigma0`	starting value for the spatial hyperparamter of CAR prior.
`psi0`	starting value for the spatial hyperparamter of CAR prior.
`Tau`	modifiable normal prior for the regression parameters with variance Tau$^2$.
`alpha`	modifiable prior distribution of hyperparamter psi (suggested values: 2, 1.5, 1, 0.5, 0).

`acceptb`	acceptance rate for the regression parameters beta.
`acceptga`	range of the acceptance rate for the spatial parameters gamma.
`acceptphi`	acceptance rate for the GP and ZIGP model specific dispersion parameter phi.
`acceptomega`	acceptance rate for the ZIP and ZIGP model specific extra proportion omega.
`acceptr`	acceptance rate for the NB model specific scale parameter r.
`acceptpsi`	acceptance rate for the spatial hyperparameter psi.
`beta`	are the parameter estimates for the regression parameters beta.
`gamma`	are the parameter estimates for the spatial parameters gamma.
`invsigsq`	are the parameter estimates for the inverse spatial hyperparameter sigma$^2$.
`psi`	are the parameter estimates for the spatial hyperparameter psi.
`phi`	are the parameter estimates for the GP and ZIGP model specific dispersion parameter phi.
`omega`	are the parameter estimates for the ZIP and ZIGP model specific extra proportion omega.
`r`	are the parameter estimates for the NB model specific scale parameter r.
`Coefficients`	are the number of parameter estimates.
`t.i`	exposure vector.

Gschloessl, Susanne (2006). Hierarchical Bayesian spatial regression models with applications to non-life insurance. Dissertation, Centre of Mathematical Sciences, Munich University of Technology, Chair in Mathematical Statistics, Munich University of Technology, Boltzmannstr. 3, D-85748 Garching near Munich.

Masterthesis: Schabenberger, Holger (2009). Spatial count regression models with applications to health insurance data. ("http://www-m4.ma.tum.de/Diplarb/").

Czado, C., Erhardt, V., Min, A., Wagner, S. (2007). Zero-inflated generalized Poisson models with regression effects on the mean, dispersion and zero-inflation level applied to patent outsourcing rates. Statistical Modelling 7 (2), 125-153.

R-package ZIGP for fitting GP, ZIP, ZIGP regression models using MLE.

data(sim.Yin)
data(sim.fm.X)
data(sim.region)
data(sim.gmat)
data(sim.nmat)
# true parameters for generating this data:
# beta.true = c(-1, 0.4, 1.5)
# gamma.true = vector of spatial effects according to the CAR model with mean 0, psi = 3 and sigma = 1
# range of gamma.true = c(-0.851, 0.8405)

# run all examples with higher number of iterations if you want to approximate the true parameters
# properly

poi <- est.sc(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region,
model="Poi", sim.gmat, sim.nmat, totalit=10)

# posterior means not considering a burn-in or thinning of iterations
apply(poi$beta,1,mean)
apply(poi$gamma,1,mean)


# Compare Poisson to different model classes
nb <- est.sc(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, model="NB", sim.gmat, sim.nmat, totalit=10)

gp <- est.sc(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, model="GP", sim.gmat, sim.nmat, totalit=10)

zip <- est.sc(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, model="ZIP", sim.gmat, sim.nmat, totalit=10)

zigp <- est.sc(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, model="ZIGP", sim.gmat, sim.nmat, totalit=10)

DIC.poi <- DIC(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, poi)
DIC.nb <- DIC(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, nb)
DIC.gp <- DIC(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, gp)
DIC.zip <- DIC(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, zip)
DIC.zigp <- DIC(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, zigp)

ll.poi <- LogLike(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, poi)
ll.nb <- LogLike(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, nb)
ll.gp <- LogLike(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, gp)
ll.zip <- LogLike(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, zip)
ll.zigp <- LogLike(sim.Yin, ~1+sim.fm.X[,1]+sim.fm.X[,2], sim.region, zigp)

Vuong.poi.nb <- Vuongtest(ll.poi, ll.nb, alpha = 0.05, p = DIC.poi$p.D, q = DIC.nb$p.D, correction = TRUE)
Vuong.poi.gp <- Vuongtest(ll.poi, ll.gp, alpha = 0.05, p = DIC.poi$p.D, q = DIC.gp$p.D, correction = TRUE)
Vuong.poi.zip <- Vuongtest(ll.poi, ll.zip, alpha = 0.05, p = DIC.poi$p.D, q = DIC.zip$p.D, correction = TRUE)
Vuong.poi.zigp <- Vuongtest(ll.poi, ll.zigp, alpha = 0.05, p = DIC.poi$p.D, q = DIC.zigp$p.D, correction = TRUE)

Clarke.poi.nb <- Clarketest(ll.poi, ll.nb, alpha = 0.05, p = DIC.poi$p.D, q = DIC.nb$p.D, correction = TRUE)
Clarke.poi.gp <- Clarketest(ll.poi, ll.gp, alpha = 0.05, p = DIC.poi$p.D, q = DIC.gp$p.D, correction = TRUE)
Clarke.poi.zip <- Clarketest(ll.poi, ll.zip, alpha = 0.05, p = DIC.poi$p.D, q = DIC.zip$p.D, correction = TRUE)
Clarke.poi.zigp <- Clarketest(ll.poi, ll.zigp, alpha = 0.05, p = DIC.poi$p.D, q = DIC.zigp$p.D, correction = TRUE)

acceptb/(i+1)  0.9166667 0.8333333 0.9166667
acceptga1/i acceptga2/(i+1)   0.6666667 0.9166667
acceptpsi/(i+1)   0.5
[1] -0.7997907  0.3782731  1.2812049
  [1] -0.131338699 -0.201964374 -0.190623189  0.193766129  0.598721261
  [6]  0.088630399 -0.680114078 -0.328106152  0.085668511 -0.170340613
 [11] -0.548369147 -0.831306382 -0.347765942  0.238895114  0.286159826
 [16] -0.151763809 -0.576205636 -0.069807763 -0.020328358 -0.446687604
 [21]  0.040326823 -0.240899747 -0.566873459 -0.072296893 -0.055072450
 [26] -0.101161540  0.103604771  0.134979584 -0.021134658 -0.227337712
 [31] -0.788359607 -0.368180622 -0.096313862  0.067101270 -0.059353784
 [36] -0.316669553 -0.227926857  0.001125308 -0.096710830 -0.422851558
 [41] -0.087885753 -0.627144994 -0.313013140  0.020420532  0.112484706
 [46] -0.184709362 -0.015330038  0.003083964 -0.276705916 -0.054352711
 [51] -0.397984265 -0.083195102 -0.535856014 -0.556866286 -0.314617129
 [56] -0.422710511 -0.173376067  0.009778570 -0.574588246 -0.051204989
 [61] -0.419713588 -0.569189486 -0.280708286 -0.264141339 -0.294156086
 [66] -0.358970798 -0.182056639 -0.033749899 -0.408795662 -0.313527683
 [71] -0.636864567 -0.547122105  0.041916377 -0.138096341  0.319944113
 [76] -0.066043354 -0.447881748  0.105171761 -0.147067586 -0.305077660
 [81] -0.355808548 -0.425615445 -0.420293498  0.198870398  0.123559042
 [86] -0.134589136 -0.036342396  0.185293673  0.301146507  0.106380064
 [91] -0.121079638 -0.309751734 -0.502860512 -0.479382267 -0.001944016
 [96] -0.449077472  0.139164112  0.225666035  0.162114934  0.204067341
acceptb/(i+1)   0.8333333 0.9166667 0.9166667
acceptga1/i acceptga2/(i+1)   0.75 0.9166667
acceptr/(i+1)  0.8333333
acceptpsi/(i+1)   0.5
acceptb/(i+1)  0.9166667 0.9166667 0.9166667
acceptga1/i acceptga2/(i+1)   0.75 0.9166667
acceptphi/(i+1)  0.9166667
acceptpsi/(i+1)   0.4166667
acceptb/(i+1)   0.9166667 0.9166667 0.8333333
acceptga1/i acceptga2/(i+1)   0.25 0.9166667
acceptomega/(i+1)  0.8333333
acceptpsi/(i+1)   0.3333333
acceptb/(i+1)   0.8333333 0.9166667 0.8333333
acceptga1/i acceptga2/(i+1)   0.4166667 0.9166667
acceptphi/(i+1)   0.9166667
acceptomega/(i+1)  0.9166667
acceptpsi/(i+1)   0.8333333
DIC             9011.955
mean deviance   8774.775
p.D             237.18
DIC             9061.449
mean deviance   8892.97
p.D             168.4787
DIC             9122.461
mean deviance   8864.102
p.D             258.3596
DIC             9020.271
mean deviance   8776.152
p.D             244.1195
DIC             9009.874
mean deviance   8779.648
p.D             230.2265
Favour model 1   0
No decision      0
Favour model 2   1
Favour model 1   1
No decision      0
Favour model 2   0
Favour model 1   0.75
No decision      0.25
Favour model 2   0
Favour model 1   0
No decision      0.1666667
Favour model 2   0.8333333
Favour model 1   0
No decision      0
Favour model 2   1
Favour model 1   1
No decision      0
Favour model 2   0
Favour model 1   0.8333333
No decision      0.1666667
Favour model 2   0
Favour model 1   0.08333333
No decision      0.5
Favour model 2   0.4166667