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
|
#
# effectfun.R
#
# $Revision: 1.25 $ $Date: 2020/01/23 05:36:51 $
#
effectfun <- local({
okclasses <- c("ppm", "kppm", "lppm", "dppm", "rppm", "profilepl")
effectfun <- function(model, covname, ..., se.fit=FALSE, nvalues=256) {
if(!inherits(model, okclasses))
stop(paste("First argument 'model' should be a fitted model of class",
commasep(sQuote(okclasses), " or ")),
call.=FALSE)
orig.model <- model
model <- as.ppm(model)
dotargs <- list(...)
#' determine names of covariates involved
intern.names <-
if(is.marked.ppm(model)) c("x", "y", "marks") else c("x", "y")
needed.names <- variablesinformula(rhs.of.formula(formula(model)))
#' check for clashes/quirks
if("lambda" %in% needed.names) {
if(is.dppm(orig.model) && (
identical.formulae(formula(model), ~offset(log(lambda))-1) ||
identical.formulae(formula(model), ~log(lambda)-1)
))
stop("effectfun is not defined for a DPP model with fixed intensity",
call.=FALSE)
intensityname <- setdiff(c("Lambda", "intensity"), needed.names)[1]
} else intensityname <- "lambda"
## validate the relevant covariate
if(missing(covname) || is.null(covname)) {
mc <- model.covariates(model)
if(length(mc) == 1) covname <- mc else stop("covname must be provided")
}
if(!is.character(covname))
stop("covname should be a character string", call.=FALSE)
if(length(covname) != 1L)
stop("covname should be a single character string", call.=FALSE)
# check that fixed values for all other covariates are provided
given.covs <- names(dotargs)
if(any(uhoh <- !(needed.names %in% c(given.covs, covname)))) {
nuh <- sum(uhoh)
stop(paste(ngettext(nuh,
"A value for the covariate",
"Values for the covariates"),
commasep(dQuote(needed.names[uhoh])),
"must be provided (as",
ngettext(nuh, "an argument", "arguments"),
"to effectfun)"))
}
#' establish type and range of covariate values
check.1.integer(nvalues)
stopifnot(nvalues >= 128)
N0 <- nvalues
if(covname == "x") {
covtype <- "real"
W <- as.owin(data.ppm(model))
Zr <- W$xrange
Zvals <- seq(from=Zr[1L], to=Zr[2L], length.out=N0)
} else if(covname == "y") {
covtype <- "real"
W <- as.owin(data.ppm(model))
Zr <- W$yrange
Zvals <- seq(from=Zr[1L], to=Zr[2L], length.out=N0)
} else if(covname == "marks") {
covtype <- "factor"
Zvals <- levels(marks(data.ppm(model)))
} else {
# covariate is external
if(is.data.frame(covdf <- model$covariates) && (covname %in% names(covdf))) {
Z <- covdf[,covname]
covtype <- typeof(Z)
if(covtype == "double")
covtype <- "real"
switch(covtype,
real={
Zr <- range(Z)
Zvals <- seq(from=Zr[1L], to=Zr[2L], length.out=N0)
},
integer={
Zr <- range(Z)
Zvals <- seq(from=Zr[1L], to=Zr[2L], by=ceiling((diff(Zr)+1)/N0))
},
factor={
Zvals <- levels(Z)
},
logical={
Zvals <- c(FALSE, TRUE)
},
stop(paste("Cannot handle covariate of type", dQuote(covtype)))
)
} else {
Z <- getdataobjects(covname,
environment(formula(model)),
model$covariates)[[1L]]
if(is.null(Z))
stop(paste("Cannot find covariate", sQuote(covname)),
call.=FALSE)
# convert to image
if(!is.im(Z))
Z <- as.im(Z, W=as.owin(model))
covtype <- Z$type
switch(covtype,
real={
Zr <- summary(Z)$range
Zvals <- seq(from=Zr[1L], to=Zr[2L], length.out=N0)
},
factor={
Zvals <- levels(Z)
},
logical={
Zvals <- c(FALSE, TRUE)
},
stop(paste("Cannot handle covariate of type", dQuote(covtype)))
)
}
}
# set up data frames of fake data for predict method
# First set up default, constant value for each covariate
N <- length(Zvals)
fakeloc <- resolve.defaults(dotargs,
list(x=0, y=0))[c("x","y")]
if(is.marked.ppm(model)) {
if("marks" %in% names(dotargs)) {
fakeloc$marks <- dotargs$marks
dotargs <- dotargs[names(dotargs) != "marks"]
} else {
lev <- levels(marks(data.ppm(model)))
fakeloc$marks <- lev[1L]
}
}
fakeloc <- lapply(fakeloc, padout, N=N)
fakecov <- lapply(dotargs, padout, N=N)
# Overwrite value for covariate of interest
if(covname %in% intern.names)
fakeloc[[covname]] <- Zvals
else fakecov[[covname]] <- Zvals
# convert to data frame
fakeloc <- do.call(data.frame, fakeloc)
fakecov <- if(length(fakecov) > 0) do.call(data.frame, fakecov) else NULL
#
# Now predict
pred <- predict(orig.model, locations=fakeloc, covariates=fakecov, se=se.fit)
if(!se.fit) lambda <- pred else {
lambda <- pred$estimate
se <- pred$se
sedf <- data.frame(se =se,
hi = lambda + 2 * se,
lo = lambda - 2 * se)
}
#
dfin <- if(!is.null(fakecov)) cbind(fakeloc, fakecov) else fakeloc
dfin <- dfin[covname]
dflam <- data.frame(lambda=lambda)
names(dflam) <- intensityname
df <- cbind(dfin, dflam)
#
if(covtype != "real") {
result <- df
if(se.fit) result <- cbind(result, sedf)
} else {
bc <- paren(covname)
result <- fv(df, argu=covname,
ylab=substitute(lambda(X),
list(X=as.name(covname),
lambda=as.name(intensityname))),
labl=c(covname,
paste("hat(%s)", bc)),
valu=intensityname, alim=Zr,
desc=c(paste("value of covariate", covname),
"fitted intensity"),
fname=intensityname)
if(se.fit) {
result <- bind.fv(result, sedf,
labl=c(paste("se[%s]", bc),
paste("%s[hi]", bc),
paste("%s[lo]", bc)),
desc=c("standard error of fitted trend",
"upper limit of pointwise 95%% CI for trend",
"lower limit of pointwise 95%% CI for trend"))
fvnames(result, ".") <- c(intensityname, "hi", "lo")
fvnames(result, ".s") <- c("hi", "lo")
formula(result) <- paste(". ~ ", covname)
}
}
return(result)
}
padout <- function(x,N) { rep.int(x[1L],N) }
effectfun
})
|
