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
|
#
# psstA.R
#
# Pseudoscore residual for unnormalised F (area-interaction)
#
# $Revision: 1.8 $ $Date: 2022/01/04 05:30:06 $
#
################################################################################
#
psstA <- function(object, r=NULL, breaks=NULL, ...,
model=NULL,
trend=~1, interaction=Poisson(),
rbord=reach(interaction), ppmcorrection="border",
correction="all",
truecoef=NULL, hi.res=NULL,
nr=spatstat.options("psstA.nr"),
ngrid=spatstat.options("psstA.ngrid")) {
if(is.ppm(object))
fit <- object
else if(is.ppp(object) || is.quad(object)) {
# convert to quadscheme
if(is.ppp(object))
object <- quadscheme(object, ...)
# fit model
if(!is.null(model))
fit <- update(model, Q=object, forcefit=TRUE)
else if(ppmcorrection == "border")
fit <- ppm(object,
trend=trend, interaction=interaction,
rbord=rbord, forcefit=TRUE)
else
fit <- ppm(object,
trend=trend, interaction=interaction,
correction=ppmcorrection, forcefit=TRUE)
} else
stop("object should be a fitted point process model or a point pattern")
rfixed <- !is.null(r) || !is.null(breaks)
# Extract data and quadrature points
Q <- quad.ppm(fit, drop=FALSE)
X <- data.ppm(fit)
U <- union.quad(Q)
Z <- is.data(Q) # indicator data/dummy
# E <- equalsfun.quad(Q)
# WQ <- w.quad(Q) # quadrature weights
# integrals will be restricted to quadrature points
# that were actually used in the fit
# USED <- getglmsubset(fit)
if(fit$correction == "border") {
rbord <- fit$rbord
b <- bdist.points(U)
USED <- (b > rbord)
bX <- bdist.points(X)
USEDX <- (bX > rbord)
} else {
USED <- rep.int(TRUE, U$n)
USEDX <- rep.int(TRUE, X$n)
}
# basic statistics
Win <- Window(X)
npts <- npoints(X)
areaW <- area(Win)
lambda <- npts/areaW
# determine breakpoints for r values
rmaxdefault <- rmax.rule("F", Win, lambda)
if(rfixed)
breaks <- handle.r.b.args(r, breaks, Win, rmaxdefault=rmaxdefault)
else {
# create fairly coarse 'r' values
r <- seq(0, rmaxdefault, length=nr)
breaks <- breakpts.from.r(r)
}
rvals <- breaks$r
rmax <- breaks$max
# residuals
res <- residuals(fit, type="raw", drop=FALSE,
new.coef=truecoef, quad=hi.res)
#
rescts <- with(res, "continuous")
# absolute weight for continuous integrals
wc <- -rescts
# initialise fv object
df <- data.frame(r=rvals, theo=0)
desc <- c("distance argument r", "value 0 corresponding to perfect fit")
ans <- fv(df, "r", substitute(bold(R)~Delta~V[A](r), NULL),
"theo", . ~ r,
alim=c(0, rmax), c("r","%s[theo](r)"), desc,
fname="bold(R)~Delta~V[A]")
#
# for efficiency, compute the largest value of distance transform
Dmax <- 0
for(i in 1:npts) {
Di <- distmap(X[-i])
Dimax <- summary(Di)$max
Dmax <- max(Dmax, Dimax)
}
Rmax <- min(max(rvals), Dmax * 1.1)
nontrivial <- (rvals <= Rmax)
trivialzeroes <- numeric(sum(!nontrivial))
# pseudosum
Ax <- areaLoss.grid(X, rvals[nontrivial], subset=USEDX, ngrid=ngrid)
C1 <- apply(Ax, 2, sum)
C1 <- c(C1, trivialzeroes)
# pseudocompensator
OK <- USED & !Z
Au <- areaGain.grid(U[OK], X, rvals[nontrivial], W=Win, ngrid=ngrid)
lamu <- matrix(wc[OK], nrow=nrow(Au), ncol=ncol(Au))
C2 <- apply(lamu * Au, 2, sum)
C2 <- c(C2, trivialzeroes)
# pseudoscore residual
Ctot <- C1 - C2
# tack on
ans <- bind.fv(ans,
data.frame(dat=C1,
com=C2,
res=Ctot),
c("Sigma~Delta~V[A](r)", "bold(C)~Delta~V[A](r)", "%s(r)"),
c("data pseudosum (contribution to %s)",
"model pseudocompensator (contribution to %s)",
"pseudoscore residual %s"),
"res")
#
# pseudovariance
# (skipped if called by envelope() etc)
#
if(correction == "all") {
lamX <- matrix(wc[USED & Z], nrow=nrow(Ax), ncol=ncol(Ax))
Var <- apply(lamu * Au^2, 2, sum) + apply(lamX * Ax^2, 2, sum)
Var <- c(Var, trivialzeroes)
# two-sigma limits
TwoSig <- 2 * sqrt(Var)
# tack on
ans <- bind.fv(ans,
data.frame(var=Var,
up=TwoSig,
lo=-TwoSig),
c("bold(C)^2~Delta~V[A](r)",
"%s[up](r)", "%s[lo](r)"),
c("pseudovariance of %s",
"upper 2sigma critical limit for %s",
"lower 2sigma critical limit for %s"),
"res")
fvnames(ans, ".") <- c("res", "up", "lo", "theo")
}
unitname(ans) <- unitname(fit)
#
return(ans)
}
|
