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.... <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@listas.c3sl.ufpr.br <mailto:r-br@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@listas.c3sl.ufpr.br <mailto:R-br@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@listas.c3sl.ufpr.br <mailto:R-br@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