For my project, I am trying to find a solution to a single value that will summarise the total degree of clustering in a landscape. My idea is to use the area under the curve between the Kiso value and the estimated Kpois to summarise the amount of clustering in a landscape. So using spatstat I first generate a landscape using the Matern Cluster Process.
library(spatstat)
a<-rMatClust(kappa=3,r=0.1,mu=50)
plot(a)
This produces a point pattern process which looks like this:
Then I calculate the k function for this landscape:
plot(Kest(a))
Producing something like this:
To integrate the object I desire I create a new object and then using the object select function I select the iso curve.
a1<-Kest(a)
a1$iso
This gives the values of Kiso, like this:
[1] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 5.347451e-05 5.347451e-05 1.336863e-04 2.673725e-04
[9] 3.208470e-04 3.893240e-04 4.427985e-04 5.497475e-04 5.497475e-04 6.566965e-04 7.903828e-04 9.775435e-04
[17] 1.057755e-03 1.137967e-03 1.218179e-03 1.378602e-03 1.565763e-03 1.699449e-03 1.839898e-03 2.027059e-03
[25] 2.160745e-03 2.347906e-03 2.750274e-03 3.044384e-03 3.338493e-03 3.391968e-03 3.593308e-03 3.807206e-03
[33] 4.074578e-03 4.368688e-03 4.720071e-03 5.014180e-03 5.335027e-03 5.556981e-03 5.797616e-03 5.904565e-03
[41] 6.118463e-03 6.546259e-03 6.917126e-03 7.344922e-03 7.505346e-03 7.993509e-03 8.427558e-03 8.849891e-03
[49] 9.365021e-03 9.739342e-03 1.002475e-02 1.047832e-02 1.090860e-02 1.140798e-02 1.180013e-02 1.204077e-02
[57] 1.237526e-02 1.278305e-02 1.338868e-02 1.387758e-02 1.429013e-02 1.464173e-02 1.526999e-02 1.578158e-02
[65] 1.610819e-02 1.657016e-02 1.689999e-02 1.735750e-02 1.754466e-02 1.811175e-02 1.857947e-02 1.903887e-02
[73] 1.977850e-02 2.028189e-02 2.088888e-02 2.149992e-02 2.220413e-02 2.282796e-02 2.339736e-02 2.404440e-02
[81] 2.455658e-02 2.513234e-02 2.550129e-02 2.635066e-02 2.696561e-02 2.757747e-02 2.825938e-02 2.889327e-02
[89] 2.934851e-02 3.027668e-02 3.088764e-02 3.144002e-02 3.202824e-02 3.260617e-02 3.343151e-02 3.403603e-02
[97] 3.472754e-02 3.521852e-02 3.582390e-02 3.614769e-02 3.672896e-02 3.743762e-02 3.795639e-02 3.878466e-02
[105] 3.918747e-02 3.987037e-02 4.056020e-02 4.142285e-02 4.230633e-02 4.309149e-02 4.397746e-02 4.457770e-02
[113] 4.517624e-02 4.601898e-02 4.679228e-02 4.760429e-02 4.838268e-02 4.895086e-02 4.969699e-02 5.067618e-02
[121] 5.124806e-02 5.229468e-02 5.277757e-02 5.333491e-02 5.411873e-02 5.495489e-02 5.575142e-02 5.646889e-02
[129] 5.729154e-02 5.800508e-02 5.867400e-02 5.947592e-02 6.027026e-02 6.103556e-02 6.197163e-02 6.272775e-02
[137] 6.372453e-02 6.426049e-02 6.484804e-02 6.547513e-02 6.642026e-02 6.716063e-02 6.813582e-02 6.892019e-02
[145] 6.977850e-02 7.035644e-02 7.125824e-02 7.204221e-02 7.283380e-02 7.378071e-02 7.493134e-02 7.572249e-02
[153] 7.626541e-02 7.711025e-02 7.799507e-02 7.871646e-02 7.961540e-02 8.051056e-02 8.130331e-02 8.206135e-02
[161] 8.329835e-02 8.423484e-02 8.494931e-02 8.585844e-02 8.655054e-02 8.743821e-02 8.823509e-02 8.883506e-02
[169] 8.987727e-02 9.069811e-02 9.179832e-02 9.256096e-02 9.350054e-02 9.435087e-02 9.531505e-02 9.615118e-02
[177] 9.704846e-02 9.820476e-02 9.902785e-02 9.984852e-02 1.007876e-01 1.017041e-01 1.029741e-01 1.040420e-01
[185] 1.050362e-01 1.059676e-01 1.068518e-01 1.080314e-01 1.090975e-01 1.101988e-01 1.112008e-01 1.121691e-01
[193] 1.133207e-01 1.140431e-01 1.148378e-01 1.156454e-01 1.167004e-01 1.177383e-01 1.187564e-01 1.195542e-01
[201] 1.204226e-01 1.212549e-01 1.222205e-01 1.231546e-01 1.236793e-01 1.247121e-01 1.254543e-01 1.263916e-01
[209] 1.273371e-01 1.281239e-01 1.291322e-01 1.296093e-01 1.301229e-01 1.306640e-01 1.314903e-01 1.326400e-01
[217] 1.336800e-01 1.343824e-01 1.350751e-01 1.357783e-01 1.368061e-01 1.375290e-01 1.385450e-01 1.392193e-01
[225] 1.400911e-01 1.409824e-01 1.422904e-01 1.429186e-01 1.437285e-01 1.448333e-01 1.457631e-01 1.465312e-01
[233] 1.473684e-01 1.481809e-01 1.492364e-01 1.502789e-01 1.513047e-01 1.520370e-01 1.528260e-01 1.535681e-01
[241] 1.547475e-01 1.557022e-01 1.563688e-01 1.571513e-01 1.580008e-01 1.587216e-01 1.596034e-01 1.603538e-01
[249] 1.614469e-01 1.621929e-01 1.631378e-01 1.644734e-01 1.653108e-01 1.659698e-01 1.668125e-01 1.677499e-01
[257] 1.684751e-01 1.695380e-01 1.704921e-01 1.714112e-01 1.723993e-01 1.732953e-01 1.739233e-01 1.747298e-01
[265] 1.755402e-01 1.764933e-01 1.774505e-01 1.783552e-01 1.792308e-01 1.801794e-01 1.813804e-01 1.822910e-01
[273] 1.831720e-01 1.839408e-01 1.845575e-01 1.855722e-01 1.862344e-01 1.872080e-01 1.879182e-01 1.885910e-01
[281] 1.895579e-01 1.903136e-01 1.913740e-01 1.920141e-01 1.927692e-01 1.937709e-01 1.946348e-01 1.956823e-01
[289] 1.965854e-01 1.971467e-01 1.977666e-01 1.986277e-01 1.997155e-01 2.007456e-01 2.013980e-01 2.020538e-01
[297] 2.028140e-01 2.037131e-01 2.047504e-01 2.055348e-01 2.063109e-01 2.071692e-01 2.079541e-01 2.087922e-01
[305] 2.096543e-01 2.104788e-01 2.112694e-01 2.118718e-01 2.129244e-01 2.138709e-01 2.146658e-01 2.155264e-01
[313] 2.162086e-01 2.171919e-01 2.180573e-01 2.191216e-01 2.199635e-01 2.205530e-01 2.212991e-01 2.221186e-01
[321] 2.230608e-01 2.239266e-01 2.246927e-01 2.256478e-01 2.266203e-01 2.272242e-01 2.280025e-01 2.287606e-01
[329] 2.294767e-01 2.298512e-01 2.306597e-01 2.314828e-01 2.319722e-01 2.326508e-01 2.331897e-01 2.339028e-01
[337] 2.347861e-01 2.355942e-01 2.361183e-01 2.370703e-01 2.375817e-01 2.384708e-01 2.392515e-01 2.399976e-01
[345] 2.407872e-01 2.418672e-01 2.425124e-01 2.432793e-01 2.440492e-01 2.445367e-01 2.456282e-01 2.464520e-01
[353] 2.473659e-01 2.483027e-01 2.490313e-01 2.497979e-01 2.503788e-01 2.511798e-01 2.519527e-01 2.528309e-01
[361] 2.538637e-01 2.545978e-01 2.556648e-01 2.561737e-01 2.569160e-01 2.577535e-01 2.583107e-01 2.590845e-01
[369] 2.598702e-01 2.607302e-01 2.613032e-01 2.620710e-01 2.628885e-01 2.638368e-01 2.645643e-01 2.655404e-01
[377] 2.663641e-01 2.672732e-01 2.681936e-01 2.688166e-01 2.698429e-01 2.707996e-01 2.713808e-01 2.720574e-01
[385] 2.726791e-01 2.733599e-01 2.742664e-01 2.747397e-01 2.753911e-01 2.759569e-01 2.765006e-01 2.770614e-01
[393] 2.777941e-01 2.786974e-01 2.794753e-01 2.800314e-01 2.806355e-01 2.808734e-01 2.815856e-01 2.828376e-01
[401] 2.833634e-01 2.843512e-01 2.854397e-01 2.862206e-01 2.870055e-01 2.877996e-01 2.883812e-01 2.889118e-01
[409] 2.896342e-01 2.904345e-01 2.913330e-01 2.919967e-01 2.927732e-01 2.935698e-01 2.941481e-01 2.948494e-01
[417] 2.955340e-01 2.965466e-01 2.973846e-01 2.981585e-01 2.990708e-01 2.998979e-01 3.005317e-01 3.011067e-01
[425] 3.020069e-01 3.026795e-01 3.037582e-01 3.046263e-01 3.056040e-01 3.063861e-01 3.071301e-01 3.079488e-01
[433] 3.082977e-01 3.091864e-01 3.099524e-01 3.109822e-01 3.117594e-01 3.126243e-01 3.134249e-01 3.144489e-01
[441] 3.155363e-01 3.164684e-01 3.177154e-01 3.184448e-01 3.193024e-01 3.201500e-01 3.209492e-01 3.219184e-01
[449] 3.226493e-01 3.233938e-01 3.241982e-01 3.248228e-01 3.256813e-01 3.264547e-01 3.273417e-01 3.282202e-01
[457] 3.292304e-01 3.304351e-01 3.313856e-01 3.323992e-01 3.333394e-01 3.340439e-01 3.351532e-01 3.363603e-01
[465] 3.374001e-01 3.385804e-01 3.392726e-01 3.398536e-01 3.409402e-01 3.418129e-01 3.425816e-01 3.434427e-01
[473] 3.444975e-01 3.455439e-01 3.466121e-01 3.474947e-01 3.481351e-01 3.491510e-01 3.499102e-01 3.511237e-01
[481] 3.517769e-01 3.529237e-01 3.538277e-01 3.548461e-01 3.559198e-01 3.572212e-01 3.581639e-01 3.590776e-01
[489] 3.604420e-01 3.618757e-01 3.628696e-01 3.637896e-01 3.649634e-01 3.657484e-01 3.666734e-01 3.675945e-01
[497] 3.683975e-01 3.692571e-01 3.703407e-01 3.716068e-01 3.725040e-01 3.732096e-01 3.739694e-01 3.751397e-01
[505] 3.763022e-01 3.773503e-01 3.782108e-01 3.790302e-01 3.801908e-01 3.812602e-01 3.825241e-01 3.835340e-01
[513] 3.845260e-01
I try to integrate these values between 0 and 0.25 using:
integrate(a1$iso,0,0.25)
But this returns the error "Error in match.fun(f) : 'a1$iso' is not a function, character or symbol"
How can I integrate the area under the curve (AUC) for the k function values?
Short answer:
integrate(as.function(a1), 0, 0.25)
Long answer:
The help for integrate
says that its first argument should be a function in the R language, representing the function that you want to integrate.
Your printout shows that a1$iso
is a vector of numbers. These are the values of the K function (isotropic correction estimate) K(r) at a sequence of values of the distance r. But the r values have been discarded, so there's no way to determine the area under the curve from a1$iso
alone.
The object a1
represents a function, except that it is really just a table of values of the function, and belongs to the class "fv"
(as you can read in the help for Kest
).
Good news --- there is a method to convert this object to a function in the R language. See the help for as.function.fv
. So you can just type
integrate(as.function(a1), 0, 0.25)
Here 0
and 0.25
are the limits on the horizontal axis, between which you want the integral to be computed. You will have to decide these limits. To do it automatically you could use something like with(a1, range(r))
Warning: the term AUC is usually reserved for the area under the Receiver Operating Characteristic (ROC) curve.