On Mon, 24 Apr 2023 12:37:39 -0400
Christopher Ryan via R-sig-Geo <r-sig-geo@r-project.org> wrote:
> Hello. I have a multi-type point pattern, a ppp object in spatstat,
> with two types, cases and controls. About 2600 point altogether.
> Overall, about 2% are cases. All in a polygonal window. I'm
> interested in the spatially varying relative risk of being a case
> rather than a control. I'm aware of the relrisk() command in
> spatstat, and also the risk() command in the sparr package, both of
> which take a ppp object as an argument. I've done it with both. They
> yield essentially similar pixel images, except that the range of the
> relative risk, as shown in the colored legend, is very different in
> the two plots.
>
> plot(relrisk(thin.ems.ppp, at = "pixels", weights = NULL, varcov =
> NULL, relative=TRUE, adjust=1, edge=TRUE, diggle=TRUE, se=FALSE,
> casecontrol=TRUE, case = "case"))
>
> yields relative risks ranging from 0 to about 0.16
>
> In contrast,
>
> risk(thin.ems.ppp, log = FALSE, verbose = TRUE, doplot = TRUE)
>
> yields an image with relative risks ranging from 0 to about 40.
>
> I've read the documentation for the packages, but perhaps I am still
> misunderstanding what each package means by "relative risk." Can
> anyone comment?
I have asked the maintainer of the sparr package (Tilman Davies) about
this issue, and after a bit of to-ing and fro-ing this is what has
turned up.
The basic reason for the discrepancy is the somewhat perplexing fact
that ‘relrisk’ computes the ratio of intensities, but ‘risk’ computes
the ratio of densities. I'm afraid that I can give you no insight in
respect of the intuitive interpretation of the two forms of "relative
risk".
Another, less fundamental, problem is that the default bandwidths
differ between the two functions. To get consistent results, these
bandwidths need to be set. (It is not entirely clear from the help,
to go about doing so.)
The following example, provided by Tilman, will provide you with some
guidance:
library(sparr)
data(pbc)
cas <- split(pbc)[[1]]
con <- split(pbc)[[2]]
# Using risk() with the "bandwidth" (the standard deviation of the
# kernel density estimator) set equal to 3 (h0 = 3).
r1 <- risk(cas,con,h0=3,log=FALSE)$rr
# "Raw" calculation, imitating what risk() does, i.e.
# using a ratio of densities.
fhat <- density.ppp(cas,sigma=3)
ghat <- density.ppp(con,sigma=3)
r2 <- (fhat/integral(fhat))/(ghat/integral(ghat))
# Using relrisk() with the "bandwidth" set equal to 3 (sigma = 3).
# Note the difference in argument names and the fact that we have to
# specify "relative=TRUE" (see the help for relrisk as well as having
# to specify which mark corresponds to the "control"
r3 <- relrisk(pbc,sigma=3,relative=TRUE,control=2)
# "Raw" calculation, imitating what relrisk() does, i.e.
# using a ratio of intensities.
r4 <- fhat/ghat
# Now view the results.
plot(anylist(risk=r1,risk.raw=r2,relrisk=r3,
relrisk.raw=r4),nrows=2,main="")
All is in harmony! OMMMMMMMMMM!
Hope this helps.
cheers,
Rolf Turner
--
Honorary Research Fellow
Department of Statistics
University of Auckland
Stats. Dep't. (secretaries) phone:
+64-9-373-7599 ext. 89622
Home phone: +64-9-480-4619
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
