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

Wagner Wolff wwolff em usp.br
Quarta Maio 10 13:55:29 BRT 2017 

Olá Elias obrigado pela ajuda!

Estou considerando só o efeito espacial (aleatório)? Pensei  que estivesse
incluso o fixo (covaráveis) também. Como acrescento os dois na predição?
Talvez isso resolva meu problema.

Quanto ao stk.val eu criei ele para comparar os dados observados com os
preditos, como explicado em
*http://www.math.sciences.univ-nantes.fr/~lavancie/slides_GT/INLA_report2012.pdf
<http://www.math.sciences.univ-nantes.fr/~lavancie/slides_GT/INLA_report2012.pdf>
*na página 25.

Outra questão é a seguinte. Os dados são do tipo proporção, nesse caso é
adequado usar a beta no inla(family=...)?

Abraço

2017-05-10 13:33 GMT-03:00 Elias T. Krainski via R-br <
r-br em listas.c3sl.ufpr.br>:

> 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='Den
> sity')
>
> 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
>
>
> _______________________________________________
> R-br mailing listR-br em listas.c3sl.ufpr.brhttps://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.
>
>
>
> _______________________________________________
> 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.
>



-- 
*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/20170510/ced1c626/attachment.html>

Mais detalhes sobre a lista de discussão R-br