geoR
. Leitura de dados, formato geodata e análises exploratórias.geoR
com o comando data(package="geoR")e efetuar análises descritivas atentando as características a serem observadas comentadas em aula (distribuição e assimetria, dados atípicos, tendências com coordenadas e covariáveis, possível padrão espacial)
install.packages("geoR", dep=T)
)grf()
sparseM
)Atividades:
Atividades:
sp
grf()
require(geoR) ml <- likfit(s100, ini=c(1, 0.3)) gr <- expand.grid(seq(0,1,len=50), seq(0,1, l=50)) args(output.control) ## definindo simulacoes nos resultados (output) OC <- output.control(n.pred=1000, simulations.pred=T) kc <- krige.conv(s100, loc=gr, krige=krige.control(obj.m=ml), out=OC) ## vendo o que tem nos resultados names(kc) str(kc) ## os simulacoes ficam armazenadas aqui dim(kc$simulations) ## calculando (predizendo) FUNCIONAIS ## FUNCIONAL 1: mapa de probabilidade do atributo estar acima de 1,8 p1.8 <- apply(kc$simulations, 1, function(x) mean(x>1.8)) length(p1.8) image(kc, val=p1.8) ## adicionando legenda args(legend.krige) legend.krige(x.leg=c(1.05, 1.1), val=p1.8, y.leg=c(0.2, 0.8), vert=T) ## mudando os limites da image para incluir a legenda image(kc, val=p1.8, xlim=c(0, 1.2)) legend.krige(x.leg=c(1.05, 1.1), val=p1.8, y.leg=c(0.2, 0.8), vert=T) ## Outro funcional: proposção da área com valores acima de 1,8 A1.8 <- apply(kc$simulations, 2, function(x) mean(x>1.8)) length(A1.8) hist(A1.8, prob=T) lines(density(A1.8)) rug(A1.8) summary(A1.8) ## outro funcional : distribuição dos máximos sobre a área MAX <- apply(kc$simulations, 2, function(x) max(x)) length(MAX) hist(MAX, prob=T) lines(density(MAX)) rug(MAX) ## probabilidade do MAX está acima de 4 mean(MAX > 4) ## idem para minimo MIN <- apply(kc$simulations, 2, function(x) min(x)) summary(MIN) hist(MIN, prob=T) lines(density(MIN)) rug(MIN) ## ## um outro funcional diferente do anterior seria um mapa de maximos POR PIXEL MAX.map <- apply(kc$simulations, 1, function(x) max(x)) length(MAX.map) image(kc, val=MAX.map)
## require(geoR) sata(s100) args(krige.bayes) args(model.control) MC <- model.control() args(prior.control() ## ## definindo prior para beta e fixando os valores dos demais parâmetros ## PC <- prior.control(beta="flat", sigmasq.pr="fixed", sigmasq=1, phi.prior="fixed", phi=0.3) ## definindo grid de predição gr <- as.matrix(expand.grid(seq(0,1,len=50), seq(0,1,len=50))) ## obtendo posterioris e preditivas s100.kb <- krige.bayes(s100, loc=gr, model=MC, prior=PC) ## inspecionando o output names(s100.kb) names(s100.kb$posterior) ## vendo os resultados da posterioris s100.kb$posterior ## e as predicoes na preditiva... names(s100.kb$predict) s100.kb$predict$mean[1:10] s100.kb$predict$var[1:10] s100.kb$predict$dist image(s100.kb, col=terrain.colors(21)) image(s100.kb, val=sqrt(s100.kb$predict$var), col=terrain.colors(21)) ## mapa de um funcional: probabilidade de estar acima de 2.0 ## calculando as probabilidades p2.0 <- pnorm(2.0, s100.kb$predictive$mean, sqrt(s100.kb$predictive$variance), low=F) ## e colocando no mapa image(s100.kb, val=p2.0, xlim=c(0, 1.2)) args(legend.krige) legend.krige(x.leg=c(1.05, 1.15), y.leg=c(0.2, 0.8), p2.0, vert=T, off=0.3) image(s100.kb, val=p2.0, xlim=c(0, 1.2), col=terrain.colors(5)) ## ## priori para beta e sigmasq ## args(prior.control) PCsig <- prior.control(beta="flat", sigmasq.pr="rec", phi.prior="fixed", phi=0.3) s100.kb.sig <- krige.bayes(s100, loc=gr, model=MC, prior=PCsig) names(s100.kb.sig) s100.kb.sig$posterior names(s100.kb.sig$predictive) ## pribalilidade na t (tem que corrigir o comando abaixo!!!!! p2.0 <- pt(2.0, s100.kb$predictive$mean, sqrt(s100.kb$predictive$variance), low=F) s100.kb.sig$predictive$dist ## ## prioris em beta, sigmasq e phi ## args(prior.control) PCphi <- prior.control(beta="flat", sigmasq.pr="rec", phi.prior="rec", phi.disc=seq(0, 1.5, by=0.1)) ## a chamada seria como abaixo (mas pode demorar muito para fazer predicao ## em um grid muito fino s100.kb.phi <- krige.bayes(s100, loc=gr, model=MC, prior=PCphi) ## definindo um grid mais "grosseiro" para teste gr <- as.matrix(expand.grid(seq(0,1,len=30), seq(0,1,len=30))) s100.kb.phi <- krige.bayes(s100, loc=gr, model=MC, prior=PCphi) names(s100.kb.names) names(s100.kb.phi) s100.kb.sig$posterior s100.kb.phi$posterior names(s100.kb.phi$posterior) s100.kb.phi$posterior$beta s100.kb.phi$posterior$sigmasq s100.kb.phi$posterior$phi names(s100.kb.phi$posterior) s100.kb.phi$posterior$sample ## visualizando as posterioris (marginais) ## beta|y hist(s100.kb.phi$posterior$sample[,1], prob=T) density(s100.kb.phi$posterior$sample[,1]) lines(density(s100.kb.phi$posterior$sample[,1])) rug(s100.kb.phi$posterior$sample[,1]) ## sigmasq|y hist(s100.kb.phi$posterior$sample[,2], prob=T, main=expression(sigma^2)) lines(density(s100.kb.phi$posterior$sample[,2])) rug(s100.kb.phi$posterior$sample[,2]) ## phi|y barplot(table(s100.kb.phi$posterior$sample[,3])) ## grafico "automatico" da geoiR com priori e posteriori plot(s100.kb.phi) ## experimentando com diferentes prioris ## note que aqui nao vamos fazer predição!!! PCphi <- prior.control(beta="flat", sigmasq.pr="rec", phi.prior="unif", phi.disc=seq(0, 1.5, by=0.05)) s100.kb.phi <- krige.bayes(s100, model=MC, prior=PCphi) plot(s100.kb.phi) ## PCphi <- prior.control(beta="flat", sigmasq.pr="rec", phi.prior="squar", phi.disc=seq(0, 1.5, by=0.05)) s100.kb.phi <- krige.bayes(s100, model=MC, prior=PCphi) plot(s100.kb.phi) ## priori com amis pontos na discretizacao PCphi <- prior.control(beta="flat", sigmasq.pr="rec", phi.prior="unif", phi.disc=seq(0, 1.5, by=0.05)) s100.kb.phi <- krige.bayes(s100, model=MC, prior=PCphi) plot(s100.kb.phi) ## especificando uma priori particular do usuário args(dgamma) curve(dgamma(x, 2, sc=0.05), from=0, to=1.5) curve(dgamma(x, 2, sc=0.1), from=0, to=1.5) curve(dgamma(x, 2, sc=0.15), from=0, to=1.5) ## discretizando PRIORphi <- dgamma(seq(0, 1.5, by=0.05), 2, sc=0.15) ## .. e garantindo que soma 1 na discreta PRIORphi <- PRIORphi/sum(PRIORphi) PCphi <- prior.control(beta="flat", sigmasq.pr="rec", phi.prior=PRIORphi, phi.disc=seq(0, 1.5, by=0.05)) s100.kb.phi <- krige.bayes(s100, model=MC, prior=PCphi) plot(s100.kb.phi)
profund020a.txt
profund020a.txt
data(parana)
) segundo o GLGM Poisson de média , , e offset pelos dados no arquivo (mesma ordem das coordenadas das estações):
(a)
Poisson homogêneo, (b)
Poisson não homogêneo, ©
Log-Cox Gaussiano. No ultimo caso verificar o efeito dos parâmetros do modelo.(a)
com marcas e pontos independentes, (b)
marcas e pontos dependentes