这个问题有点陈旧，但我忍不住将@ryanm 评论（我刚刚注意到）当作一个有趣的挑战。正如我在上面的评论中提到的，问题在于 plotK 函数如何操纵设备。在 plotK 函数中对（不必要的？）代码进行一些修剪可以解决问题：
---
title: "Untitled"
output: html_document
---
```{r}
library(stpp)
data(fmd)
data(northcumbria)
FMD<-as.3dpoints(fmd[,1]/1000,fmd[,2]/1000,fmd[,3])
Northcumbria=northcumbria/1000
# estimation of the temporal intensity
Mt<-density(FMD[,3],n=1000)
mut<-Mt$y[findInterval(FMD[,3],Mt$x)]*dim(FMD)[1]
# estimation of the spatial intensity
h<-mse2d(as.points(FMD[,1:2]), Northcumbria, nsmse=50, range=4)
h<-h$h[which.min(h$mse)]
Ms<-kernel2d(as.points(FMD[,1:2]), Northcumbria, h, nx=5000, ny=5000)
atx<-findInterval(x=FMD[,1],vec=Ms$x)
aty<-findInterval(x=FMD[,2],vec=Ms$y)
mhat<-NULL
for(i in 1:length(atx)) mhat<-c(mhat,Ms$z[atx[i],aty[i]])
# estimation of the STIK function
u <- seq(0,10,by=1)
v <- seq(0,15,by=1)
stik1 <- STIKhat(xyt=FMD, s.region=northcumbria/1000,t.region=c(1,200),
lambda=mhat*mut/dim(FMD)[1], dist=u, times=v, infectious=TRUE)
```
```{r,echo=FALSE}
plotK <- function (K, n = 15, L = FALSE, type = "contour", legend = TRUE,
which = NULL, main = NULL, ...)
{
old.par <- par(no.readonly = TRUE)
on.exit(par(old.par))
correc = c("none", "isotropic", "border", "modified.border",
"translate")
correc2 = K$correction
id <- match(correc2, correc, nomatch = NA)
if ((is.null(which) && length(id) > 1) || any(is.na(match(which,
correc, nomatch = NA)))) {
mess <- paste("Please specify the argument 'which', among:",
paste(dQuote(correc2), collapse = ", "))
stop(mess, call. = FALSE)
}
if (isTRUE(K$infectious))
which = "isotropic"
if (is.matrix(K$Khat)) {
if (is.null(which))
which = correc2
else {
if (!(is.null(which)) && which != correc2) {
mess <- paste("Argument 'which' should be", paste(dQuote(correc2),
collapse = ", "))
stop(mess, call. = FALSE)
}
}
}
if (!is.matrix(K$Khat)) {
id <- match(which, correc2, nomatch = NA)
if (is.na(id)) {
mess <- paste("Please specify the argument 'which', among:",
paste(dQuote(correc2), collapse = ", "))
stop(mess, call. = FALSE)
}
else K$Khat = K$Khat[[id]]
}
if (!is.null(main)) {
titl = main
subtitl = ""
if (isTRUE(L))
k <- K$Khat - K$Ktheo
else k <- K$Khat
}
else {
if (isTRUE(L)) {
k <- K$Khat - K$Ktheo
subtitl <- paste("edge correction method: ", which,
sep = "")
if (isTRUE(K$infectious))
titl <- expression(hat(K)[ST] * group("(", list(u,
v), ")") - pi * u^2 * v)
else titl <- expression(hat(K)[ST] * group("(", list(u,
v), ")") - 2 * pi * u^2 * v)
}
else {
k <- K$Khat
titl = expression(hat(K)[ST] * group("(", list(u,
v), ")"))
subtitl <- paste("edge correction method: ", which,
sep = "")
}
}
typeplot = c("contour", "image", "persp")
id <- match(type, typeplot, nomatch = NA)
if (any(nbg <- is.na(id))) {
mess <- paste("unrecognised plot type:", paste(dQuote(type[nbg]),
collapse = ", "))
stop(mess, call. = FALSE)
}
if ((length(id) != 1) || is.na(id))
stop("Please specify one type among \"contour\", \"image\" and \"persp\" ")
typeplot = rep(0, 3)
typeplot[id] = 1
colo <- colorRampPalette(c("red", "white", "blue"))
M <- max(abs(range(k)))
M <- pretty(c(-M, M), n = n)
n <- length(M)
COL <- colo(n)
if (typeplot[3] == 1) {
mask <- matrix(0, ncol = length(K$times), nrow = length(K$dist))
for (i in 1:length(K$dist)) {
for (j in 1:length(K$times)) {
mask[i, j] <- COL[findInterval(x = k[i, j], vec = M)]
}
}
COL <- mask[1:(length(K$dist) - 1), 1:(length(K$times) -
1)]
if (isTRUE(legend)) {
par(cex.lab = 2, cex.axis = 1.5, font = 2, lwd = 1,
mar = c(0, 0, 3, 0))
par(fig = c(0, 0.825, 0, 1))
persp(x = K$dist, y = K$times, z = k, xlab = "u",
ylab = "v", zlab = "", expand = 1, col = COL,
...)
title(titl, cex.main = 1.5, sub = subtitl, outer = TRUE,
line = -1)
par(fig = c(0.825, 1, 0, 1))
mini <- findInterval(x = min(k, na.rm = TRUE), vec = M)
maxi <- findInterval(x = max(k, na.rm = TRUE), vec = M)
legend("right", fill = colo(n)[maxi:mini], legend = M[maxi:mini],
horiz = F, bty = "n")
}
else {
par(cex.lab = 2, cex.axis = 1.5, font = 2, lwd = 1)
persp(x = K$dist, y = K$times, z = k, xlab = "u",
ylab = "v", zlab = "", expand = 1, col = COL,
...)
title(titl, cex.main = 1.5, sub = subtitl)
}
}
if (typeplot[1] == 1) {
if (isTRUE(legend)) {
par(cex.lab = 1.5, cex.axis = 1.5, font = 2, plt = c(0,
1, 0, 1), lwd = 1, mar = c(0.5, 0.5, 2.5, 0.5),
las = 1)
par(fig = c(0.1, 0.825, 0.1, 1))
contour(K$dist, K$times, k, labcex = 1.5, levels = M,
drawlabels = F, col = colo(n), zlim = range(M),
axes = F)
box(lwd = 2)
at <- axTicks(1)
axis(1, at = at[1:length(at)], labels = at[1:length(at)])
at <- axTicks(2)
axis(2, at = at[1:length(at)], labels = at[1:length(at)])
title(titl, cex.main = 1.5, sub = subtitl, outer = TRUE,
line = -1)
par(fig = c(0, 1, 0.1, 1))
mini <- findInterval(x = min(k, na.rm = TRUE), vec = M)
maxi <- findInterval(x = max(k, na.rm = TRUE), vec = M)
legend("right", fill = colo(n)[maxi:mini], legend = M[maxi:mini],
horiz = F, bty = "n")
}
else {
par(cex.lab = 2, cex.axis = 1.5, font = 2, lwd = 2,
las = 1)
contour(K$dist, K$times, k, labcex = 1.5, levels = M,
drawlabels = T, col = colo(n), zlim = range(M),
axes = F)
box(lwd = 2)
at <- axTicks(1)
axis(1, at = at[1:length(at)], labels = at[1:length(at)])
at <- axTicks(2)
axis(2, at = at[1:length(at)], labels = at[1:length(at)])
title(titl, cex.main = 1.5, sub = subtitl)
}
}
if (typeplot[2] == 1) {
if (isTRUE(legend)) {
par(cex.lab = 1.5, cex.axis = 1.5, font = 2, lwd = 1,
plt = c(0, 1, 0, 1), mar = c(0.5, 0.5, 2.5, 0.5),
las = 1)
par(fig = c(0.1, 0.825, 0.1, 1))
image(K$dist, K$times, k, col = colo(n), zlim = range(M),
axes = F, xlab = "", ylab = "")
box(lwd = 2)
at <- axTicks(1)
axis(1, at = at[1:length(at)], labels = at[1:length(at)])
at <- axTicks(2)
axis(2, at = at[1:length(at)], labels = at[1:length(at)])
title(titl, cex.main = 1.5, sub = subtitl, outer = TRUE,
line = -1)
par(fig = c(0, 1, 0.1, 1))
mini <- findInterval(x = min(k, na.rm = TRUE), vec = M)
maxi <- findInterval(x = max(k, na.rm = TRUE), vec = M)
legend("right", fill = colo(n)[maxi:mini], legend = M[maxi:mini],
horiz = F, bty = "n")
}
else {
par(cex.lab = 2, cex.axis = 1.5, font = 2, lwd = 2,
las = 1)
image(K$dist, K$times, k, col = colo(n), zlim = range(M),
axes = F, xlab = "", ylab = "")
box(lwd = 2)
at <- axTicks(1)
axis(1, at = at[1:length(at)], labels = at[1:length(at)])
at <- axTicks(2)
axis(2, at = at[1:length(at)], labels = at[1:length(at)])
title(titl, cex.main = 1.5, sub = subtitl)
}
}
par(old.par)
}
```
```{r}
plotK(stik1)
```
如果您经常使用 stpp 包，可能值得向维护人员发送一封电子邮件，说明为什么需要弄乱设备。