Hi All
Hope someone can help.
I'm trying to run some spatial models on a binomial and a poisson
dataset, and so geoRglm seemed like a good choice.
I'm also running the same models in SAS with proc Glimmix. Since SAS
uses a quasiliklihood I get no AIC values for model ranking, and I'm
able to calculate them from R readily enough given the number of
parameters etc.
My problem is this: The mcmc model requires a prior beta in order to
estimate parameters. The parameter estimates I get never seem to vary
far from my initial estimates, and the resulting AIC is heavily
influenced by this. So if I start all betas at 0 I get one answer,
and if I seed it with the parameter estimates from a SAS model I get
a reflection in the parameter estimates of what I start it with, but
VERY different model rankings. I have posted my attempt at code, and
model results using SAS estimates for priors, and 0s as estimates for
priors below.
Can anyone shed light on how to get parameter estimates to converge
on a value less influenced by my choice of priors, and whether or not
my ranking method is appropriate?
Thanks for your help.
Ken
#Transet length, habitat model and single track roads using SAS priors
WMmcmccontrl4 <- mcmc.control(S.scale = .3, thin=100, phi.start=0.15,
phi.scale = 0.1)
WMliktrend4 <- trend.spatial(~TRAN_LNTH + GRSPSTBG + STAll, WMgeo)
WMmcmcmodel4 <- list(cov.pars = c(.1067,.15), beta = c
(-2.739,0.0003,0.0275,-0.0310), family="binomial",
link='logit',trend= WMliktrend4)
WMmcmctest4 <- glsm.mcmc(WMgeo, coords=WMgeo$coords, model =
WMmcmcmodel4, mcmc.input=WMmcmccontrl4)
#have a look at convergence of the chain etc.
WMmcmctest4.c <- create.mcmc.coda(WMmcmctest4$simulations[45,],
mcmc.input=WMmcmccontrl4)
par(mfrow = c(1,2))
plot(WMmcmctest4.c, density=FALSE, ask=FALSE, auto.layout =FALSE)
autocorr.plot(WMmcmctest4.c, ask=FALSE, auto.layout=FALSE)
#fit parameters
WMmcmcprep4 <- prepare.likfit.glsm(WMmcmctest4)
WMmcmclikfit4 <- likfit.glsm(WMmcmcprep4, cov.model="spherical",
trend= WMliktrend4, ini.phi=.1, fix.nugget.rel = TRUE,
limits=pars.limits(phi=c(0.01,10)))
WMmcmclikfit4
summary(WMmcmclikfit4)
# calculate AIC
WMmcmclikfit4.AIC <- -2*WMmcmclikfit4$loglik + 2* WMmcmclikfit4$npars
WMmcmclikfit4.AIC
> WMmcmclikfit4
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi
"-2.6830" " 0.0003" " 0.0268" "-0.0312" " 0.0925" " 0.1000"
likfit.glsm : maximised log-likelihood = 2.713
> WMmcmclikfit4.AIC
[1] 6.573307
> summary(WMmcmclikfit4)
Summary of the maximum likelihood parameter estimation
-----------------------------------
Family = binomial , Link = logit
Parameters of the mean component (trend):
beta0 beta1 beta2 beta3
"-2.6830134" " 0.0003121" " 0.0267731" "-0.0312180"
Parameters of the spatial component:
correlation function: spherical
(estimated) variance parameter sigmasq (partial sill) = 0.0925
(estimated) cor. fct. parameter phi (range parameter) = 0.1
(fixed) relative nugget = 0
Maximised Likelihood:
log.L n.params
"2.713" "6"
Call:
likfit.glsm(mcmc.obj = WMmcmcprep4, trend = WMliktrend4, cov.model =
"spherical",
ini.phi = 0.1, fix.nugget.rel = TRUE, limits = pars.limits(phi =
c(0.01,
10)))
>
Using 0 for prior parameter estimates
> WMmcmclikfit4
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi
"-0.1995" " 0.0000" " 0.0028" "-0.0028" " 0.1098" " 0.1000"
likfit.glsm : maximised log-likelihood = 10.48
> summary(WMmcmclikfit4)
Summary of the maximum likelihood parameter estimation
-----------------------------------
Family = binomial , Link = logit
Parameters of the mean component (trend):
beta0 beta1 beta2 beta3
"-1.995e-01" " 1.432e-05" " 2.838e-03" "-2.793e-03"
Parameters of the spatial component:
correlation function: spherical
(estimated) variance parameter sigmasq (partial sill) = 0.1098
(estimated) cor. fct. parameter phi (range parameter) = 0.1
(fixed) relative nugget = 0
Maximised Likelihood:
log.L n.params
"10.48" "6"
Call:
likfit.glsm(mcmc.obj = WMmcmcprep4, trend = WMliktrend4, cov.model =
"spherical",
ini.phi = 0.1, fix.nugget.rel = TRUE, limits = pars.limits(phi =
c(0.01,
10)))
> # calculate AIC
> WMmcmclikfit4.AIC <- -2*WMmcmclikfit4$loglik + 2* WMmcmclikfit4$npars
> WMmcmclikfit4.AIC
[1] -8.969521
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
