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

Quarta Maio 10 13:33:05 BRT 2017

Olá Wagner,

No seu cenário de predição você está considerando apenas o efeito 
espacial. Isso só fará sentido num cenário sem covariáveis ou quando 
nenhuma for importante.

Bwt, não entendi porque você repete exatamente o mesmo cenário em 
'stk.est' e 'stk.val'... para mim é redundante.

Att.

Elias.

On 09/05/2017 17:18, Wagner Wolff via R-br wrote:
> 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 <tel:+55%2019%2098238-5582>
> http://orcid.org/0000-0003-3426-308X
> https://github.com/wwolff7
> http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1
>
>
> _______________________________________________
> R-br mailing list
> R-br em listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20170510/4917ff49/attachment.html>