1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
|
# Copyright (c) 2007-2008 Markus Reder and Roger Bivand
lm.morantest.exact <- function(model, listw, zero.policy = NULL,
alternative = "greater", spChk=NULL, resfun=weighted.residuals,
zero.tol=1.0e-7, Omega=NULL, save.M=NULL, save.U=NULL, useTP=FALSE,
truncErr=1e-6, zeroTreat=0.1)
{
if (!inherits(listw, "listw"))
stop(paste(deparse(substitute(listw)), "is not a listw object"))
if (!inherits(model, "lm"))
stop(paste(deparse(substitute(model)), "not an lm object"))
if (is.null(zero.policy))
zero.policy <- get("zeroPolicy", envir = .spdepOptions)
stopifnot(is.logical(zero.policy))
N <- length(listw$neighbours)
u <- resfun(model)
if (N != length(u))
stop("objects of different length")
if (is.null(spChk)) spChk <- get.spChkOption()
if (spChk && !chkIDs(u, listw))
stop("Check of data and weights ID integrity failed")
if (!(alternative %in% c("greater", "less", "two.sided")))
stop("alternative must be one of: \"greater\", \"less\", or \"two.sided\"")
u <- as.vector(u)
listw.U <- listw2U(listw)
S0 <- sum(unlist(listw.U$weights))
lu <- lag.listw(listw.U, u, zero.policy = zero.policy)
Nnn <- N
if (zero.policy) Nnn <- length(which(card(listw$neighbours) > 0))
I <- (Nnn/S0) * (crossprod(u, lu) / crossprod(u))
I_save <- I
if (!isTRUE(all.equal((Nnn/S0), 1))) I <- I * (S0/Nnn)
p <- model$rank
p1 <- 1:p
nacoefs <- which(is.na(coefficients(model)))
XtXinv <- chol2inv(model$qr$qr[p1, p1, drop = FALSE])
X <- model.matrix(terms(model), model.frame(model))
# fixed after looking at TOWN dummy in Boston data
if (length(nacoefs) > 0L) X <- X[,-nacoefs]
if (!is.null(wts <- weights(model))) {
X <- sqrt(diag(wts)) %*% X
}
M <- diag(N) - X %*% tcrossprod(XtXinv, X)
U <- listw2mat(listw.U)
if (is.null(Omega)) {
MVM <- M %*% U %*% M
MVM <- 0.5 * (t(MVM) + MVM)
evalue <- eigen(MVM, only.values=TRUE)$values
idxpos <- which(abs(evalue) < zero.tol)
if (length(idxpos) != p)
warning("number of zero eigenvalues greater than number of variables")
idxpos <- idxpos[1] - 1
if (idxpos < 1) {
warning("first eigenvalue index zero")
gamma <- c()
} else gamma <- evalue[1:idxpos]
gamma <- c(gamma, evalue[(idxpos+1+p):N])
res <- exactMoran(I, gamma, alternative=alternative, type="Global",
useTP=useTP, truncErr=truncErr, zeroTreat=zeroTreat)
} else {
if (dim(Omega)[1] != N) stop("Omega of different size than data")
res <- exactMoranAlt(I, M, U, Omega, N, alternative=alternative,
type="Alternative", useTP=useTP, truncErr=truncErr,
zeroTreat=zeroTreat)
}
data.name <- paste("\nmodel:", paste(strwrap(gsub("[[:space:]]+", " ",
paste(deparse(model$call), sep="", collapse=""))), collapse="\n"),
"\nweights: ", deparse(substitute(listw)), "\n", sep="")
res$estimate <- c(I_save)
res$data.name <- data.name
res$df <- (N-p)
if (!is.null(save.M)) res$M <- M
if (!is.null(save.U)) res$U <- U
return(res)
}
# function contributed by Michael Tiefelsdorf 2008
# implements his 2000 book eq. 6.7, p.69
truncPoint <- function(SpecI, truncErr=1e-6, zeroTreat=0.1){
m <- length(SpecI)
absSpecI <- abs(SpecI)
absSpecI[absSpecI < zeroTreat] <- zeroTreat
TU <- (truncErr*pi*m/2)^(-2/m) * prod(absSpecI^(-1/m))
# Break product up to reduce rounding errors and over- or under-flows
return(TU)
}
exactMoran <- function(I, gamma, alternative="greater", type="Global", np2=NULL, useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) {
if (!(alternative %in% c("greater", "less", "two.sided")))
stop("alternative must be one of: \"greater\", \"less\", or \"two.sided\"")
if (type == "Global") {
SpecI <- gamma - c(I)
integrand <- function(x) {
sin(0.5 * colSums(atan(SpecI %*% t(x)))) /
(x * apply((1 + SpecI^2 %*% t(x^2))^(1/4), 2, prod))}
} else if (type == "Alternative") {
if (useTP) SpecI <- gamma
integrand <- function(x) {sin(0.5 * colSums(atan((gamma) %*% t(x)))) /
(x * apply((1+(gamma)^2 %*% t(x^2))^(1/4), 2, prod))}
} else if (type == "Local") {
if (is.null(np2)) stop("Local requires np2")
min <- gamma[1]
max <- gamma[2]
if (useTP) SpecI <- c(min, rep(0,np2), max) - I
integrand <- function(x) {
sin(0.5*(atan((min-I)*x)+(np2)*atan((-I)*x) +
atan((max-I)*x)))/(x*((1+(min-I)^2*x^2) *
(1+I^2*x^2)^(np2) * (1+(max-I)^2*x^2))^(1/4))}
}
if (useTP) upper <- truncPoint(SpecI, truncErr=truncErr,
zeroTreat=zeroTreat)
else upper <- Inf
II <- integrate(integrand, lower=0, upper=upper)$value
# FIXME II > pi/2
if (II > pi/2 && type == "Local") {
tau <- gamma
df <- np2 + 2
if (length(tau) == 2L) tau <- c(tau[1], rep(0, df-2), tau[2])
E.I <- sum(tau)/df
tau <- tau - E.I
V.I <- (2*sum(tau^2)) / (df*(df+2))
sd.ex <- (I - E.I) / sqrt(V.I)
warning("Normal approximation SD substituted", call.=FALSE)
oType <- "N"
} else {
sd.ex <- qnorm(0.5-II/pi)
oType <- "E"
}
if (alternative == "two.sided") p.v <- 2 * pnorm(abs(sd.ex),
lower.tail=FALSE)
else if (alternative == "greater")
p.v <- pnorm(sd.ex, lower.tail=FALSE)
else p.v <- pnorm(sd.ex)
if (!is.finite(p.v) || p.v < 0 || p.v > 1)
warning("Out-of-range p-value: reconsider test arguments")
statistic <- sd.ex
attr(statistic, "names") <- "Exact standard deviate"
p.value <- p.v
estimate <- c(I)
attr(estimate, "names") <- "Observed Moran I"
method <- paste(type, "Moran I statistic with exact p-value")
res <- list(statistic = statistic, p.value = p.value,
estimate = estimate, method = method,
alternative = alternative, gamma=gamma, oType=oType)
class(res) <- "moranex"
return(res)
}
exactMoranAlt <- function(I, M, U, Omega, n, alternative="greater",
type="Alternative", useTP=FALSE, truncErr=1e-6, zeroTreat=0.1) {
Omega <- chol(Omega)
A <- Omega %*% M %*% (U - diag(I, n)) %*% M %*% t(Omega)
gamma <- sort(eigen(A)$values)
obj <- exactMoran(I, gamma, alternative=alternative,
type=type, useTP=useTP, truncErr=truncErr, zeroTreat=zeroTreat)
obj
}
print.moranex <- function(x, ...) {
class(x) <- c("htest", "moranex")
print(x, ...)
invisible(x)
}
H1_moments <- function(M, U, Omega, n) {
B <- Omega %*% M %*% t(Omega)
eigen <- eigen(B)
lambda <- eigen$values
P <- eigen$vectors
A <- Omega %*% M %*% U %*% M %*% t(Omega)
H <- t(P) %*% A %*% P
hh <- diag(H)
integrand <- function(x) apply((1+2*lambda %*% t(x))^(-0.5),2,prod) *
colSums(hh/(1+2*lambda %*% t(x)))
mu <- integrate(integrand,lower=0, upper=Inf)$value
integrand2 <- function(x) {
res=0
for (i in 1:n){
for(j in 1:n){
res=res+(H[i,i]*H[j,j]+2*H[i,j]^2) /
((1+2*lambda[i]*x)*(1+2*lambda[j]*x))*x
}
}
apply((1+2*lambda %*% t(x))^(-0.5),2,prod)*res
}
mu2 <- integrate(integrand2,lower=0, upper=Inf)$value
res<-list(Ew=mu,Var=mu2-mu^2)
res
}
moranExpect_H1 <- function(listw, rho, select=FALSE){
if (!(select[1]))
select=1:length(listw$neighbours)
n <- length(select)
V <- listw2mat(listw)[select,select]
M <- diag(n)-matrix(rep(1/n,n*n),nrow=n)
WOm <- invIrW(listw, rho=rho)[select,select]
B <- WOm %*% M %*% t(WOm)
eigen <- eigen(B)
lambda <- eigen$values
P <- eigen$vectors
A <- WOm %*% M %*% (0.5*(V + t(V))) %*% M %*% t(WOm)
H <- t(P) %*% A %*% P
hh <- diag(H)
integrand <- function(x) apply((1+2*lambda %*% t(x))^(-0.5), 2, prod) *
colSums(hh/(1+2*lambda %*% t(x)))
mu <- integrate(integrand, lower=0, upper=Inf)$value
mu
}
moranExpect_rho_H1 <- function(listw, I0, select=FALSE){
if (!(select[1]))
select=1:length(listw$neighbours)
V <- listw2mat(listw)[select,select]
eigenvalues<-as.double(eigen(V)$values)
min=min(eigenvalues)
max=max(eigenvalues)
f <- function(x) abs(I0-moranExpect_H1(listw, x, select=select))
rop <- optimise(f, c(1/min,1/max))$minimum
rop
}
|
