[R-br] Intervalo de predição INLA
Elias T. Krainski
eliaskrainski em yahoo.com.br
Quarta Maio 10 15:27:56 BRT 2017
Wagner,
Você não colocou covariáveis na 'stk.pred' (veja o comentário no topo da
pag. 24 do link que vc enviou).
Não é necessário fazer uma stack para validacao se não há dados fora
daqueles usados para estimar o modelo (para ser realmente uma validação)
pois os preditos para 'stk.est' e 'stk.pred' serão a mesma coisa dado
que o scenario é o mesmo. Note que na pagina 25 do link foi usado outro
scenario para validação, com 367 locais.
Sobre o suporte da resposta, a verossimilhanca beta é uma boa opção.
Elias
On 10/05/2017 13:55, Wagner Wolff wrote:
> 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/%7Elavancie/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 <mailto: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='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
>> <http://orcid.org/0000-0003-3426-308X>
>> https://github.com/wwolff7
>> http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1
>> <http://buscatextual.cnpq.br/buscatextual/visualizacv.do?id=K4463141A1>
>>
>>
>> _______________________________________________
>> R-br mailing list
>> R-br em listas.c3sl.ufpr.br <mailto:R-br em listas.c3sl.ufpr.br>
>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-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 <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 <mailto:R-br em listas.c3sl.ufpr.br>
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-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
> <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
> <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
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20170510/f91f3a0e/attachment.html>
Mais detalhes sobre a lista de discussão R-br