[R-br] Intervalo de predição INLA

[Wagner Wolff]

Olá pessoal estou fazendo uma modelagem geoestatística pelo INLA, mas estou
com dúvidas quanto às estimativas para os intervalos de credibilidade
maiores, p.ex. 95%, para esta situação os valores estimados fogem do campo
amostral que é de 0,3 a 0,7. Alguém sabe onde posso configurar para que as
estimativas fiquem nesse intervalo. Segue o código:

## Criando domain
IEBdomain <- inla.nonconvex.hull(as.matrix(dados[,1:2]), -0.03, -0.05,
resolution=c(100,100))

## Crando mesh
IEBmesh <- inla.mesh.2d(boundary=IEBdomain, max.edge=c(35,35), cutoff=35,
offset=c(-0.5, -0.5))
plot(IEBmesh, asp=1, main='')

## spde matern 0.5 = exponetial
IEBspde <- inla.spde2.matern(mesh=IEBmesh,alpha=2)

mesh.index <- inla.spde.make.index(name = "i",
                                  n.spde = IEBspde$n.spde)

## Matriz projetora estimativa
A.est <- inla.spde.make.A(IEBmesh, loc=as.matrix(dados[,1:2]))

## Matriz de covariaveis selecionadas pelo AIC, estatistica frequentista
covars <- dados[,c(1:4,6:23)]

stk.est <- inla.stack(data=list(y=dados$IEB_ANO), A=list(A.est,1),
tag="est",
                      effects=list(c(mesh.index,list(Intercept=1)),
                                     list(covars)))

stk.val <- inla.stack(data=list(y=NA), A=list(A.est,1), tag="est",
                      effects=list(c(mesh.index,list(Intercept=1)),
                                     list(covars)))
## Matriz projetora predicao
A.pred = inla.spde.make.A(IEBmesh)
stk.pred = inla.stack(data = list(y = NA),
                      A = list(A.pred),tag = "pred",
                      effects=list(c(mesh.index,list(Intercept=1))))

str(stk.pred)
stk.all <- inla.stack(stk.est, stk.val,stk.pred)

## Testar qual variável tem menor DIC
names(covars)
f.IEB <- y ~ -1 + Intercept + Dens.dren + f(i, model=IEBspde)
names(inla.models()$likelihood)
r.IEB <-inla(f.IEB,family="beta",
control.compute=list(dic=TRUE),quantiles=c(0.025,0.1,0.5,
0.975),
              data=inla.stack.data(stk.all,spde=IEBspde),
control.predictor=list(A=inla.stack.A(stk.all),compute=TRUE))

names(r.IEB)
r.IEB$dic$dic
r.IEB$summary.fixed
r.IEB$summary.hyper[1,]
r.IEB$summary.hyper[-1,]

result <- inla.spde2.result(r.IEB, "i", IEBspde)
names(result)
str(r.IEB$marginals.hyperpar)

## Posterior mean
inla.emarginal(function(x) x, result$marginals.variance.nominal[[1]])
inla.emarginal(function(x) x, result$marginals.range.nominal[[1]])

## Quantis
inla.qmarginal(c(0.025,0.5,0.975), result$marginals.variance.nominal[[1]])
inla.qmarginal(c(0.025,0.5,0.975), result$marginals.range.nominal[[1]])

par(mfrow=c(2,3), mar=c(3,3.5,0,0), mgp=c(1.5, .5, 0), las=0)

plot(r.IEB$marginals.fix[[1]], type='l', xlab=expression(beta[0]),
ylab='Density')
plot(r.IEB$marginals.fix[[2]], type='l', xlab=expression(beta[1]),ylab=
'Density')
plot(r.IEB$marginals.hy[[1]], type='l', xlab=expression(phi),ylab='Density')

plot.default(inla.tmarginal(function(x) 1/exp(x), r.IEB$marginals.hy[[3]]),
type='l',
             xlab=expression(kappa), ylab='Density')
plot.default(result$marginals.variance.nominal[[1]], type='l',
xlab=expression(sigma[x]^2), ylab='Density')
plot.default(result$marginals.range.nominal[[1]], type='l', xlab='Practical
range',
             ylab='Density')

index.pred <- inla.stack.index(stk.all, "pred")$data

names(r.IEB$summary.linear.predictor)

linpred.mean <- r.IEB$summary.linear.predictor[index.pred,"mean"]
linpred.2.5 <- r.IEB$summary.linear.predictor[index.pred,"0.025quant"]
linpred.10 <- r.IEB$summary.linear.predictor[index.pred,"0.1quant"]
linpred.50 <- r.IEB$summary.linear.predictor[index.pred,"0.5quant"]
linpred.97.5 <- r.IEB$summary.linear.predictor[index.pred,"0.975quant"]

(nxy <- round(c(diff(c(200,800)), diff(c(6700,7200)))))
proj <- inla.mesh.projector(IEBmesh, xlim=c(200,800), ylim=c(6700,7200),
dims=nxy)

lp.mean.grid <- inla.mesh.project(proj, linpred.mean)
lp.2.5.grid <- inla.mesh.project(proj, linpred.2.5)
lp.10.grid <- inla.mesh.project(proj, linpred.10)
lp.50.grid <- inla.mesh.project(proj, linpred.50)
lp.97.5.grid <- inla.mesh.project(proj, linpred.97.5)

par(mfrow=c(2,3), mar=c(3,3.5,0,0), mgp=c(1.5, .5, 0), las=0)

image(lp.2.5.grid)
image(lp.10.grid)
image(lp.mean.grid)
image(lp.97.5.grid)


Abraço

-- 
*Wagner Wolff, **PhD*
"*Luiz de Queiroz**" College of Agriculture,*
University of São Paulo
Pádua Dias avenue11 | 13418-900| Piracicaba-SP| Brazil
Phone:  +55 19 982385582 <+55%2019%2098238-5582>
http://orcid.org/0000-0003-3426-308X
https://github.com/wwolff7
http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20170509/3a52aa15/attachment.html>