Adrian, Thanks for providing a better alternative for analyzing the data. You are right in that the density of people exposed to snake bites is not uniform across the study area. I've been thinking about either incorporating human population density as a covariate in the model or correcting snake bite counts by population density.
Marcelo Marcelo Tognelli <[email protected]> writes: > I have probability maps of the distribution of 4 species of venomous snakes > (raster files output from species distribution modeling software) and point > locality data with information on snake bite events (most of them without > the id of the species involved in the accident). I would like to run an > analysis to see what species correlates best with snake bite events. > My idea is to generate a kernel density raster from the point event data and > then do > some kind of spatial correlation against the species distribution maps. You don't need to smooth the point event data to obtain a correlation. The correlation between a point process and a random field is well-defined. In the 'spatstat' package, if X is a point pattern and Z is a pixel image, then the sample correlation can be computed by Xp <- pixellate(X, W=as.owin(Z)) cor(as.vector(as.matrix(Z)), as.vector(as.matrix(Xp)), use="pairwise.complete.obs") However, the interpretation of correlations in a spatially-inhomogeneous dataset is beset with problems. Probably the best way to analyse these data is to regard the snake bite events as a point process 'response' and the species distribution data as the 'covariates' that may be used to explain the response. This is similar to many standard analyses in spatial epidemiology. The very simplest statistical assumption would be that each species has a different intrinsic probability of biting (the probability that in a particular confrontation a particular snake will bite). In that case the intensity of snake bites (number of snake bites per unit area, in the survey period) would be proportional to p1 * a1 + p2 * a2 + p3 * a3 + p4 * a4 where p1, p2..p4 are the (unknown) bite probabilities of each species, and a1, a2...a4 are the (known) spatially-varying densities of abundance for the 4 species. Assuming snakes act independently of each other (!) the snake bites are a Poisson process with this intensity. Unfortunately this model cannot (yet) be fitted in 'spatstat' because it is linear and the current implementation of 'ppm' is restricted to canonical (loglinear) models. The following hack would work. objfun <- function(beta, X, Z1, Z2, Z3, Z4) { p <- exp(beta) inten <- eval.im(p[1] * Z1 + p[2]* Z2 + p[3]*Z3 + p[4] * Z4) -sum(log(inten[X])) + summary(inten)$integral } op <- optim(rep(0,4), objfun, X=bites, Z1=species1, Z2=species2, Z3=species3, Z4=species4) p <- exp(op$par) where 'bites' is your point pattern of snake bites, and species1, .., species4 are pixel images of the distributions of the 4 snakes. The result 'p' gives maximum likelihood estimates of the relative biting propensities for the 4 species (these are not normalised so only their relative values makes sense). By tweaking the call to 'eval.im' you can try different kinds of models. There's an issue here about the 'opportunity' or 'exposure' to snake bite. If the density of people (or visits by people) is not uniform across the study region then this could affect the observed distribution of snake bites. "Tomislav Hengl" <[email protected]> writes: > Surprisingly, I could not find any function in the spatstat package (or > splancs package) that > specifically derives cross-correlations between multiple point processes. The Mark Connection Function 'markconnect' does exactly this, e.g. data(lansing) plot(alltypes(lansing, markconnect)) (But this is for correlations between different points, not between points and spatial variables.) > # fit stationary marked Poisson process with different intensity for each > species: > lansing.ppm <- ppm(lansing, ~marks, Poisson()) > summary(lansing.ppm) > ...but this does not say anything about which species are most correlated > (and which are negatively correlated). Of course not - the fitted model assumes they are NOT correlated. Adrian Baddeley _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
