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/~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@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@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@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