est.sc {spatcounts} | R Documentation |
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)