require(geoR) data(Ksat) points(Ksat) ## caso nao tenha as coordenadas da borda da area...uma possivbel forma... #plot(Ksat$coords, asp=1) #bor <- locator(type="l") #bor <- cbind(bor$x, bor$y) plot(Ksat) ## buscando a transformação require(MASS) boxcox(Ksat) plot(Ksat) plot(Ksat, lam=0) ## nota: no caso da transformação logarítimica considere tb o uso ## de logtrans() logtrans(Ksat$dat~1) ## ou seja, neste caso nao precisamos somar nada aos dados! ### outros dados data(parana) plot(parana) ## inadequado: boxcox(parana) ## adequado: buscar a transformação considerando média não cte boxcox(parana, trend="1st") ## outro dado ainda data(ca20) plot(ca20) points(ca20) polygon(ca20$reg1) polygon(ca20$reg2) polygon(ca20$reg3) summary(ca20) ## inadequado: boxcox(ca20) ## adequado: buscar a transformação considerando média não cte boxcox(ca20, trend=~area+altitude) ## nao precisa transformação (lambda =1) ## Escolhando modelo para ca20 m0 <- likfit(ca20, ini=c(100, 100), nug=20) m1 <- likfit(ca20, ini=c(100, 100), nug=20, trend=~area) m2 <- likfit(ca20, ini=c(100, 100), nug=20, trend=~area+altitude) m0 m1 m2 m0$AIC m1$AIC m2$AIC m0$BIC m1$BIC m2$BIC ## checando consistencia de resultados em relacao oas valores iniciais m1 <- likfit(ca20, ini=c(100, 100), nug=20, trend=~area) m1 m1 <- likfit(ca20, ini=c(90, 100), nug=30, trend=~area) m1 ## ou ainda... formacendo varios valores #m1 <- likfit(ca20, ini=expand.grid(seq(20,120,by=20), seq(40, 200,by=40)), nug=c(0, 10, 20, 30, 40), trend=~area) #m1 ## e o efeito espacial? vale a pena incluí-lo summary(m1)