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]]


