Modelo Poisson com transformação
10
A figura abaixo apresenta os pontos amostrados.
●
●
●
●
●
●
●
●
●
●
8
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
6
●
●
●
●
●●
●
●
●
●
●
●
●
●
●●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
2
●
●
●
4
COORD[,2]
●
●
●
●
●
●
2
4
6
8
10
COORD[,1]
Figura 1: Coordenadas amostrais
Os dados foram simulados utilizando os comandos a seguir:
>
>
>
>
>
>
>
MU=5 ; sigma2=2; phi=2
#simulando os dados utilizando distribuiç~
ao Poisson
D<-as.matrix(dist(COORD,diag=T,upper=T))
Sigma<-sigma2*exp(-D/phi)
S=t(chol(as.matrix(Sigma)))%*%rnorm(nrow(COORD))
lambdaS=exp(MU+S)
Y=rpois(nrow(COORD),lambdaS)
Após a simulação a necessidade de transformação foi verificada por meio
do gráfico boxcox, e como esperado a transformação logarı́tmica foi sugerida.
Para verficar a existência de padrão espacial utilizou-se o variograma
com envelope de simulação, apresentado na figura .
Foram empregados três modelos, cujos comandos estão a seguir.
> M1<-likfit(dad,ini=c(0.7*var(Y),0.3*var(Y)))
> M2<-likfit(dad,trend = "1st",ini=c(0.7*var(Y),0.3*var(Y)))
1
−550
−700
−750
−900
−850
−800
log−Likelihood
−650
−600
95%
−2
−1
0
1
2
λ
Figura 2: Gráfico para verificação de transformação
>
>
>
>
>
>
>
>
>
M3<-likfit(dad,trend = "2nd",ini=c(0.7*var(Y),0.3*var(Y)))
L1=-2*M1$loglik;L2=-2*M2$loglik;L3=-2*M3$loglik
pchisq(L1-L2,M2$npars-M1$npars,lower.tail=F)
pchisq(L1-L3,M3$npars-M1$npars,lower.tail=F) #Ficamos com o modelo 3
M1$AIC; M2$AIC;M3$AIC
M1R<-likfit(dad,ini=c(0.7*var(Y),0.3*var(Y)), lik.method = "REML")
M2R<-likfit(dad,trend = "1st",ini=c(0.7*var(Y),0.3*var(Y)), lik.method = "REML")
M3R<-likfit(dad,trend = "2nd",ini=c(0.7*var(Y),0.3*var(Y)), lik.method = "REML")
M1R$AIC; M2R$AIC;M3R$AIC
Dos modelos utilizados, tando pelo método da máxima verossilhança
(ML) como pelo método da máxima verossimilhança restrita o modelo selecionado foi o modelo 3, quando utilizando a semente 75. Os valores dos AIC
para cada cada um dos modelos foram 218.85, 222.03,225, utilizando ML. A
figura apresenta os mapas de krigagem do modelo selecionado para ambos
os métodos.
Além do estudo anterior foram simuladas 100 amostras de tamanho 75,
e calculados os componentes de variância pelos métodos ML e REML, como
apresentado a seguir:
2
10
10
●
●
●
● ●
8
Y Coord
6
Y Coord
6
8
●
●
●
●
●
●●
●
●
4
4
●
●
●
●
●
●
●
2
2
●
●
4
6
X Coord
8
10
2
●
●
●
●
●
●
●
●
● ●
●
●
●●
● ● ● ●●
● ● ●
●●
●
●
●
●
●
●
●
● ●
● ●
● ● ●
●
●
●● ● ●
●
●● ● ●● ● ● ●
●
●
●
● ●
●
● ● ●
●
●
●
●
●
●●
●
2
4
data
6
●
●
4
6
X Coord
●
● ●
●●
●
●
●
●●
●
● ●
●
●
●
●
●● ●
● ●●
●
●●
●
●●
●●
●
● ●
●
●
●
● ●●
●
●
●●
●●●●
●
● ● ●●
● ●
●
●
●
● ●
●
●
●
●
●
●
4
6
8
data
●
●
●
2
●
●
Density
0.00 0.05 0.10 0.15 0.20 0.25
8
2
●
●
●
8
10
2
4
6
8
data
Figura 3: Análise exploratória para os dados transformados
>
>
>
>
>
+
+
+
+
+
+
+
+
+
+
+
+
+
+
require(geoR)
NSIM <- 100
sigmasqM=sigmasqMR=phiM=phiMR=tauM=tauMR=numeric(NSIM)
MU=5 ; sigma2=2; phi=2
for (i in 1:NSIM){
COORD<-cbind(sample(seq(1,10,length=100),75),sample(seq(1,10,length=100),75))
D<-as.matrix(dist(COORD,diag=T,upper=T))
Sigma<-sigma2*exp(-D/phi)
S=t(chol(as.matrix(Sigma)))%*%rnorm(nrow(COORD))
lambdaS=exp(MU+S)
Y=rpois(nrow(COORD),lambdaS)
#usando a transformaç~
ao log
dad<-as.geodata(cbind(COORD,log(Y+1)),coords.col=c(1,2), data.col=3)
VAR.Y <- var(log(Y+1))
M3<-likfit(dad,ini=c(0.9*VAR.Y,max(D)/15), nugget=0.1*VAR.Y)
M3R<-likfit(dad,ini=c(0.9*VAR.Y,max(D)/15), nugget=0.1*VAR.Y, lik.method = "REML")
sigmasqM[i]=M3$sigmasq
sigmasqMR[i]=M3R$sigmasq
phiM[i]=M3$phi
3
4
3
2
semivariance
●
●
●
●
●
●
●
●
●
1
●
●
●
0
●
0
1
2
3
4
5
distance
Figura 4: Variograma com envelope simulado
+ phiMR[i]=M3R$phi
+ }
Após a simulação foram gerados os gráficos boxplot para os valores estimados de φ e σ 2 , figura .
Modelo Normal
Os dados foram simulados utilizando as mesmas coordenadas do modelo de
Poisson, utilizando os comandos a seguir:
>
>
>
>
>
>
>
MU=5 ; sigma2=4; phi=2 ;tau=1.5
D<-as.matrix(dist(COORD,diag=T,upper=T))
Sigma<-sigma2*exp(-D/phi)
SigmaE<-tau*diag(nrow(COORD))
E=t(chol(as.matrix(SigmaE)))%*%rnorm(nrow(COORD))
S=t(chol(as.matrix(Sigma)))%*%rnorm(nrow(COORD))
Y= MU+ S + E
4
Para verficar a existência de padrão espacial utilizou-se o variograma
com envelope de simulação, apresentado na figura .
Foram empregados três modelos, cujos comandos estão a seguir.
>
>
>
>
>
>
>
>
>
>
>
M1<-likfit(dad,ini=c(0.7*var(Y),0.3*var(Y)))
M2<-likfit(dad,trend = "1st",ini=c(0.7*var(Y),0.3*var(Y)))
M3<-likfit(dad,trend = "2nd",ini=c(0.7*var(Y),0.3*var(Y)))
L1=-2*M1$loglik;L2=-2*M2$loglik;L3=-2*M3$loglik
pchisq(L1-L2,M2$npars-M1$npars,lower.tail=F)
pchisq(L1-L3,M3$npars-M1$npars,lower.tail=F) #Ficamos com o modelo 3
M1$AIC; M2$AIC;M3$AIC
M1R<-likfit(dad,ini=c(0.7*var(Y),0.3*var(Y)), lik.method = "REML")
M2R<-likfit(dad,trend = "1st",ini=c(0.7*var(Y),0.3*var(Y)), lik.method = "REML")
M3R<-likfit(dad,trend = "2nd",ini=c(0.7*var(Y),0.3*var(Y)), lik.method = "REML")
M1R$AIC; M2R$AIC;M3R$AIC
Dos modelos utilizados, tando pelo método da máxima verossilhança
(ML) como pelo método da máxima verossimilhança restrita o modelo selecionado foi o modelo 3, quando utilizando a semente 75. Os valores dos AIC
para cada cada um dos modelos foram 320.06, 319.84,320.07, utilizando ML.
A figura apresenta os mapas de krigagem do modelo selecionado para
ambos os métodos.
Além do estudo anterior foram simuladas 100 amostras de tamanho 75,
e calculados os componentes de variância pelos métodos ML e REML, como
apresentado a seguir:
>
>
>
>
>
+
+
+
+
+
+
+
+
+
+
+
+
+
+
require(geoR)
NSIM <- 100
sigmasqM=sigmasqMR=phiM=phiMR=tauM=tauMR=numeric(NSIM)
MU=5 ; sigma2=4; phi=2 ;tau=1.5
for (i in 1:NSIM){
COORD<-cbind(sample(seq(0,10,length=100),75),sample(seq(0,10,length=100),75))
D<-as.matrix(dist(COORD,diag=T,upper=T))
Sigma<-sigma2*exp(-D/phi)
SigmaE<-tau*diag(nrow(COORD))
E=t(chol(as.matrix(SigmaE)))%*%rnorm(nrow(COORD))
S=t(chol(as.matrix(Sigma)))%*%rnorm(nrow(COORD))
#MU = as.vector(cbind(1, COORD) %*% BETA)
Y= MU+ S + E
VAR.Y <- var(Y)
dad<-as.geodata(cbind(COORD,Y),coords.col=c(1,2), data.col=3)
M3<-likfit(dad,ini=c(0.9*VAR.Y,max(D)/15), nugget=0.1*VAR.Y)
M3R<-likfit(dad,ini=c(0.9*VAR.Y,max(D)/15), nugget=0.1*VAR.Y, lik.method = "REML")
sigmasqM[i]=M3$sigmasq
sigmasqMR[i]=M3R$sigmasq
5
+
+
+
+
+
phiM[i]=M3$phi
phiMR[i]=M3R$phi
tauM[i]=M3$tau
tauMR[i]=M3R$tau
}
Após a simulação foram gerados os gráficos boxplot para os valores estimados de φ, σ 2 e τ , figura .
6
10
Método ML
●
●
●
●
●
8
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
4
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
2
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
6
●
●
●
●
Y Coord
●
●
●
●
●
●
●
0
●
0
2
4
6
8
10
X Coord
(a)
10
Método REML
●
●
●
●
●
8
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
4
●
●
●
●
●
●
●
●●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
2
●
●
●●
●
●
●
●
●
6
●
●
●
Y Coord
●
●
●
●
●
●
●
0
●
0
2
4
6
8
10
X Coord
(b)
7
Figura 5: Mapa de krigagem utilizando os métodos ML e REML
●
●
80
40
100
50
120
●
60
30
●
●
20
10
40
20
●
●
●
●
●
●
●
●
●
●
●
●
●
●
0
0
●
●
●
φ
σ2
Figura 6: Gráfico boxplot
8
10
10
●
●
●
●
●
●
●
Y Coord
6
Y Coord
6
●
●
●
●
2
●
2
●
4
4
●
●
●
4
●
6
X Coord
8
10 12
2
●●
10
0
2
4
●
●
6
data
8
4
6
data
8
●
●
10
12
●
Density
0.10
●
●
●● ● ● ●
●
●
● ●
●●
●●
●
●
● ●●
●
●
●
● ●
●● ●
●●
●
●
●●
●
●
● ●
●●
●
●●
●
● ●
●
●
●
●
●
●
●●
●
●● ●
● ●
●
●
●
0.05
8
●
●
0.15
● ●
●
●
●
0
2
●
●
●●
●
●
●
data
4
6
●
●
●
●
● ● ●
●
●
● ●
●●
●●
●
●
●
●
●
●
● ●
●
● ●●●
●
●
● ● ●●
●
●
●●●
●
●
●
●
●
●●
●
●
●
●
●● ●
●
● ●
●
●
● ●
8
8
●
●
●
●
●
●
●
2
4
6
X Coord
8
0.00
●
10
−2
0
2
Figura 7: Análise exploratória
9
10 12
8
6
4
●
●
●
●
●
●
●
●
●
2
●
●
0
semivariance
●
●
0
1
2
3
4
5
distance
Figura 8: Variograma com envelope simulado
10
10
Método ML
●
●●
●
●
●
● ●
●
●
8
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
4
●
●
●
●
●
●
●
●
●
6
Y Coord
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●●
●
●
2
●●
●
●
●
●
●
●
●
●
0
●
0
2
4
6
8
10
X Coord
(a)
10
Método REML
●
●●
●
●
●
●
● ●
●
●
8
●
●
●
● ●
●
●
6
●
●
●
●
●
●
●
●
●
●
●●
4
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
Y Coord
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
2
●
● ●
●
●
●
●
●
●
●
●
●
0
●
0
2
4
6
8
10
X Coord
(b)
11
Figura 9: Mapa de krigagem utilizando os métodos ML e REML
●
●
●
3.0
400
200
●
●
●
2.5
●
2.0
300
150
●
1.5
200
100
●
●
●
0.5
50
100
1.0
●
σ2
●
●
●
●
0.0
0
0
●
●
●
φ
Figura 10: Gráfico boxplot
12
τ