diff --git a/R/GlobalEnvelope.R b/R/GlobalEnvelope.R index 9221378..9f60384 100644 --- a/R/GlobalEnvelope.R +++ b/R/GlobalEnvelope.R @@ -1,34 +1,43 @@ -GlobalEnvelope <- -function(Simulations, Alpha) { +GlobalEnvelope <- function( + Simulations, + Alpha) { verifyclass(Simulations, "fv") # Initialize Simulations <- as.data.frame(Simulations)[, -1] # Eliminate r NumberOfSimulations <- ncol(Simulations) - TargetNumber <- (1-Alpha)*NumberOfSimulations + TargetNumber <- (1 - Alpha) * NumberOfSimulations KeptSimulations <- Simulations - PresentMinValues <- apply(Simulations, 1, min, na.rm=TRUE) - PresentMaxValues <- apply(Simulations, 1, max, na.rm=TRUE) + PresentMinValues <- apply(Simulations, 1, FUN = min, na.rm = TRUE) + PresentMaxValues <- apply(Simulations, 1, FUN = max, na.rm = TRUE) # Loop until the target number of simulations is kept while(ncol(KeptSimulations) > TargetNumber) { # Remember previous min and max PreviousMinValues <- PresentMinValues PreviousMaxValues <- PresentMaxValues # Select the simulations that gave extreme values - SimulationsToDrop <- c(unlist(apply(KeptSimulations, 1, which.min)), unlist(apply(KeptSimulations, 1, which.max))) + SimulationsToDrop <- c( + unlist(apply(KeptSimulations, MARGIN = 1, FUN = which.min)), + unlist(apply(KeptSimulations, MARGIN = 1, FUN = which.max)) + ) # Drop them KeptSimulations <- KeptSimulations[, -SimulationsToDrop] # Fails if no simulations are left - if (is.null(dim(KeptSimulations))) + if (is.null(dim(KeptSimulations))) { stop("Global envelope could not be calculated. More simulations are necessary.") + } # Calculate min and max - PresentMinValues <- apply(KeptSimulations, 1, min, na.rm=TRUE) - PresentMaxValues <- apply(KeptSimulations, 1, max, na.rm=TRUE) + PresentMinValues <- apply(KeptSimulations, 1, FUN = min, na.rm = TRUE) + PresentMaxValues <- apply(KeptSimulations, 1, FUN = max, na.rm = TRUE) } # Interpolate because the kept number of simulations is not always the target NumberOfKeptSimulations <- ncol(KeptSimulations) - Glo <- PresentMinValues + (PreviousMinValues-PresentMinValues)/NumberOfSimulations*(TargetNumber-NumberOfKeptSimulations) - Ghi <- PresentMaxValues + (PreviousMaxValues-PresentMaxValues)/NumberOfSimulations*(TargetNumber-NumberOfKeptSimulations) + Glo <- PresentMinValues + + (PreviousMinValues - PresentMinValues) / NumberOfSimulations * + (TargetNumber - NumberOfKeptSimulations) + Ghi <- PresentMaxValues + + (PreviousMaxValues - PresentMaxValues) / NumberOfSimulations * + (TargetNumber - NumberOfKeptSimulations) return(rbind(Glo, Ghi)) } diff --git a/R/GoFtest.R b/R/GoFtest.R index 9cef133..8534d1a 100644 --- a/R/GoFtest.R +++ b/R/GoFtest.R @@ -1,9 +1,9 @@ -GoFtest <- -function(Envelope) { +GoFtest <- function(Envelope) { # Verify Envelope - if (!inherits(Envelope, "envelope")) + if (!inherits(Envelope, "envelope")) { stop("Envelope is not of class envelope") + } # Verify simulations if (is.null(attr(Envelope, "simfuns"))) { stop("Envelope does not contain simulations in its attribute simfuns") @@ -14,24 +14,31 @@ function(Envelope) { } NumberOfSimulations <- dim(SimulatedValues)[2] - AverageSimulatedValues <- apply(SimulatedValues, 1, sum)/(NumberOfSimulations-1) - rIncrements <- (r-c(0,r)[seq_along(r)])[-1] + AverageSimulatedValues <- apply(SimulatedValues, 1, sum) / (NumberOfSimulations - 1) + rIncrements <- (r - c(0,r)[seq_along(r)])[-1] # Ui calculate the statistic for a simulation Ui <- function(SimulationNumber) { - Departure <- (SimulatedValues[, SimulationNumber]-AverageSimulatedValues)[seq_along(r)-1] - WeightedDeparture <- (Departure[!is.nan(Departure)])^2*rIncrements[!is.nan(Departure)] + Departure <- (SimulatedValues[, SimulationNumber] - + AverageSimulatedValues)[seq_along(r) - 1] + WeightedDeparture <- (Departure[!is.nan(Departure)])^2 * + rIncrements[!is.nan(Departure)] return(sum(WeightedDeparture)) } # Calculate the Ui statistic for all simulations - SimulatedU <- vapply(seq_len(NumberOfSimulations), Ui, FUN.VALUE=0.0) + SimulatedU <- vapply( + seq_len(NumberOfSimulations), + FUN = Ui, + FUN.VALUE = 0 + ) # Calculate the statistic for the actual value - RecenteredValues <- (ActualValues-AverageSimulatedValues)[seq_along(r)-1] - WeightedRecenteredValues <- (RecenteredValues[!is.nan(RecenteredValues)])^2*rIncrements[!is.nan(RecenteredValues)] + RecenteredValues <- (ActualValues - AverageSimulatedValues)[seq_along(r) - 1] + WeightedRecenteredValues <- (RecenteredValues[!is.nan(RecenteredValues)])^2 * + rIncrements[!is.nan(RecenteredValues)] ActualU <- sum(WeightedRecenteredValues) # Return the rank - return(mean(ActualU1) - { + invsqrtmat <- function(mat) { + if(length(mat) > 1) { e <- eigen(mat) # Eigen vectors p <- e$vectors @@ -101,11 +102,9 @@ Ktest <- function(X, r) # put into a diagonal matrix rd <- diag(d) # Resolution - return(solve(p%*%rd)) - } - else - { - return(1/sqrt(mat)) + return(solve(p %*% rd)) + } else { + return(1 / sqrt(mat)) } } #------------------------------------------------------------------------------- @@ -113,56 +112,104 @@ Ktest <- function(X, r) # r1, r2 : distances # w : width of the rectangle window # l : length of the rectangle window - covh1 <- function(r1, r2, w, l) - { - # Values of the product in the corner. Coordinates are x'=(n-xi)/r2 without normalization in r'^2r^2/w^2/l^2 - # This function must be inside covh1 because it needs to use local variables, that can not be passed as parameters because adaptIntegrate forbides it. - corner <- function(x){(foncA1(ra2*x/ra1, ra1, w, l)+foncA2(ra2*x/ra1, ra1, w, l)+foncA3(ra2*x/ra1, ra1, w, l)+foncA4(ra2*x/ra1, ra1, w, l))*(foncA3(x, ra2, w, l)+foncA4(x, ra2, w, l))} + covh1 <- function(r1, r2, w, l) { + # Values of the product in the corner. + # Coordinates are x'=(n-xi)/r2 without normalization in r'^2r^2/w^2/l^2 + # This function must be inside covh1 because it needs to use local variables, + # that can not be passed as parameters because adaptIntegrate forbides it. + corner <- function(x){ + (foncA1(ra2 * x / ra1, ra1, w, l) + + foncA2(ra2 * x / ra1, ra1, w, l) + + foncA3(ra2 * x / ra1, ra1, w, l) + + foncA4(ra2 * x / ra1, ra1, w, l)) * + (foncA3(x, ra2, w, l) + foncA4(x, ra2, w, l) + ) + } # r values must be ordered ra1 <- min(r1, r2) ra2 <- max(r1, r2) # Size ratios parameters - rapr <- ra1/ra2 - r12 <- ra1*ra1/w/l - r22 <- ra2*ra2/w/l + rapr <- ra1 / ra2 + r12 <- ra1 * ra1 / w / l + r22 <- ra2 * ra2 / w / l - # Biases, normalized by n^2/r^2 + # Biases, normalized by n^2 / r^2 b1 <- brn(ra1, w, l) b2 <- brn(ra2, w, l) # Numerical computing of the elliptic integral - int2 <- stats::integrate(integrand3, lower=0, upper=1, r1=ra1, r2=ra2) - intcorner <- cubature::adaptIntegrate(corner, lowerLimit=c(0, 0), upperLimit=c(1, 1)) + int2 <- stats::integrate( + integrand3, + lower = 0, + upper = 1, + r1 = ra1, + r2 = ra2 + ) + intcorner <- cubature::adaptIntegrate( + corner, + lowerLimit = c(0, 0), + upperLimit = c(1, 1) + ) # line 1 - covh1_ <- (w-2*r2)/w*(l-2*r2)/l*b1*b2 + covh1_ <- (w - 2 * r2) / w * (l - 2 * r2) / l * b1 * b2 # line 2 - covh1_ <- covh1_+2*(w+l-4*r2)*r2/w/l*b1*(b2-foncG(1)) + covh1_ <- covh1_ + 2 * (w + l - 4 * r2) * r2 / w / l * b1 * (b2 - foncG(1)) # line 3 - covh1_ <- covh1_+2*(w+l-4*r2)*r1/w/l*(int2$value-b2*foncG(1))+4*r22*intcorner$integral + covh1_ <- covh1_ + + 2 * (w + l - 4 * r2) * r1 / w / l * (int2$value - b2 * foncG(1)) + + 4 * r22 * intcorner$integral # multiplication by the common factor - covh1_ <- covh1_*r12*r22 + covh1_ <- covh1_ * r12 * r22 # Result return(covh1_) } - foncg <- function(x){if (x<=1) { return(acos(x)-x*sqrt(1-x*x))} else {return(0)}} - foncg01 <- function(x){acos(x)-x*sqrt(1-x*x)} - foncG <- function(x){x*acos(x)-sqrt(1-x*x)*(2+x*x)/3+2/3} + foncg <- function(x){ + if (x<=1) { + return(acos(x) - x * sqrt(1 - x * x)) + } else { + return(0) + } + } + foncg01 <- function(x){ + acos(x) - x * sqrt(1 - x * x) + } + foncG <- function(x){ + x * acos(x) - sqrt(1 - x * x) * (2 + x * x) / 3 + 2 / 3 + } # Values of h1 on different zones without normalization by r^2/w/l - brn <- function(r, w, l){4*r*(w+l)/(3*w*l)-r*r/(2*w*l)} - indic <- function(a){as.numeric(a)} - foncA1 <- function(x, r, w, l){brn(r, w, l)*indic(x[1] >= 1)*indic(x[2] >= 1)} - foncA2 <- function(x, r, w, l){(brn(r, w, l)-foncg(x[2]))*indic(x[1] >= 1)*indic(x[2] < 1)+(brn(r, w, l)-foncg(x[1]))*indic(x[2] >= 1)*indic(x[1] < 1)} - foncA3 <- function(x, r, w, l){(brn(r, w, l)-foncg(x[1])-foncg(x[2]))*indic(x[1] < 1)*indic(x[2] < 1)*indic(x[1]^2+x[2]^2 > 1)} - foncA4 <- function(x, r, w, l){(brn(r, w, l)+x[1]*x[2]-(foncg(x[1])+foncg(x[2]))/2-pi/4)*indic(x[1] < 1)*indic(x[2] < 1)*indic(x[1]^2+x[2]^2 <= 1)} + brn <- function(r, w, l){ + 4 * r * (w + l) / (3 * w * l) - r * r / (2 * w * l) + } + indic <- function(a){ + as.numeric(a) + } + foncA1 <- function(x, r, w, l){ + brn(r, w, l) * indic(x[1] >= 1) * indic(x[2] >= 1) + } + foncA2 <- function(x, r, w, l){ + (brn(r, w, l) - foncg(x[2])) * indic(x[1] >= 1) * indic(x[2] < 1) + + (brn(r, w, l) - foncg(x[1])) * indic(x[2] >= 1) * indic(x[1] < 1) + } + foncA3 <- function(x, r, w, l){ + (brn(r, w, l) - foncg(x[1]) - + foncg(x[2])) * indic(x[1] < 1) * indic(x[2] < 1) * indic(x[1]^2 + x[2]^2 > 1) + } + foncA4 <- function(x, r, w, l){ + (brn(r, w, l) + x[1] * x[2] - + (foncg(x[1]) + foncg(x[2])) / 2 - pi / 4) * indic(x[1] < 1) * + indic(x[2] < 1) * indic(x[1]^2 + x[2]^2 <= 1) + } - integrand3 <- function(x, r1, r2){foncg01(r1*x/r2)*foncg01(x)} + integrand3 <- function(x, r1, r2){ + foncg01(r1 * x / r2) * foncg01(x) + } #------------------------------------------------------------------------------- # End of utilities @@ -170,30 +217,37 @@ Ktest <- function(X, r) if (!is.rectangle(X$window)) stop("The window must be a rectangle to apply Ktest.") # value = NA if the number of points is less than 2 - if (X$n > 1) - { + if (X$n > 1) { # Width - Length - l <- X$window$xrange[2]-X$window$xrange[1] - w <- X$window$yrange[2]-X$window$yrange[1] - # espKi : expectation of K, rho unknown, calculated according to the number of points. + l <- X$window$xrange[2] - X$window$xrange[1] + w <- X$window$yrange[2] - X$window$yrange[1] + # espKi : expectation of K, rho unknown, + # calculated according to the number of points. espKi_ <- espKi(r, X$n, w, l) # Computed variance. Depends on the the observed number of points and the window. sigmaKi_ <- sigmaKi(r, X$n, w, l) # Distance matrix. - pairdist_ <- pairdist.ppp(X) # Requires a lot of RAM. Limited to 8,000 points with 32-bit versions of R. + # Requires a lot of RAM. Limited to 8,000 points with 32-bit versions of R. + pairdist_ <- pairdist.ppp(X) # Estimation of K - # NbPairs : number of pairs of points less than r apart (it's a vector, one value for each r) + # NbPairs : number of pairs of points less than r apart + # (it's a vector, one value for each r) # pairdist_ > 0 eliminates distance from a point to itself - # *1.0 is to convert values to real and avoid infinite values in NbPairs*w*l later - NbPairs <- vapply(r, function(d) sum(pairdist_ > 0 & pairdist_ < d)*1.0, FUN.VALUE=0.0) + # * 1.0 is to convert values to real and avoid infinite values + # in NbPairs * w * l later + NbPairs <- vapply( + r, + FUN = function(d) sum(pairdist_ > 0 & pairdist_ < d) * 1.0, + FUN.VALUE = 0.0 + ) # Kest gets the estimator of K, centered on the expected value. - Kest <- NbPairs*w*l/(X$n*(X$n-1))-espKi_ + Kest <- NbPairs * w * l / (X$n * (X$n - 1)) - espKi_ TestVector <- invsqrtmat(sigmaKi_) %*% Kest - return(1 - stats::pchisq(sum(TestVector*TestVector), length(r))) + return(1 - stats::pchisq(sum(TestVector * TestVector), length(r))) } else return(NA) } diff --git a/R/Smooth.wmppp.R b/R/Smooth.wmppp.R index 882ceaf..b65abb5 100644 --- a/R/Smooth.wmppp.R +++ b/R/Smooth.wmppp.R @@ -22,7 +22,7 @@ Smooth.wmppp <- function(X, fvind, distance = NULL, Quantiles = FALSE, # Reduce the point pattern to the reference type if (ReferenceType != "") { - is_ReferenceType <- spatstat.geom::marks(X)$PointType == ReferenceType + is_ReferenceType <- marks(X)$PointType == ReferenceType X <- X[is_ReferenceType] } @@ -39,7 +39,7 @@ Smooth.wmppp <- function(X, fvind, distance = NULL, Quantiles = FALSE, r_to_plot <- max(fvind$r[fvind$r<=distance]) # Weights if (Weighted) { - weights <- spatstat.geom::marks(X)$PointWeight + weights <- marks(X)$PointWeight } else { weights <- rep(1, X$n) } @@ -53,11 +53,11 @@ Smooth.wmppp <- function(X, fvind, distance = NULL, Quantiles = FALSE, if (Quantiles) { # Smooth the quantiles of the dbm # Make the quantiles the marks of X - spatstat.geom::marks(X) <- Qvalues + marks(X) <- Qvalues # Smooth() requires the top class of X to be ppp class(X) <- "ppp" # Eliminate NA's before smoothing - is_na <- is.na(spatstat.geom::marks(X)) + is_na <- is.na(marks(X)) weights <- weights[!is_na] X<- X[!is_na] Image <- Smooth.ppp(X, sigma = sigma, ..., weights = weights, adjust = Adjust, dimyx = c(Nby, Nbx)) @@ -65,11 +65,11 @@ Smooth.wmppp <- function(X, fvind, distance = NULL, Quantiles = FALSE, # Smooth the values of the dbm fvind.matrix <- as.matrix(fvind) # Extract the values. Columns 1 to 3 contain the global dbm - spatstat.geom::marks(X) <- fvind.matrix [which(fvind.matrix [, 1] == r_to_plot), -(1:3)] + marks(X) <- fvind.matrix [which(fvind.matrix [, 1] == r_to_plot), -(1:3)] # Smooth requires the top class of X to be ppp class(X) <- "ppp" # Eliminate NA's before smoothing - is_na <- is.na(spatstat.geom::marks(X)) + is_na <- is.na(marks(X)) weights <- weights[!is_na] X<- X[!is_na] Image <- Smooth.ppp(X, sigma = sigma, ..., weights = weights, adjust = Adjust, dimyx = c(Nby, Nbx)) diff --git a/R/as.Dtable.data.frame.R b/R/as.Dtable.data.frame.R index da76a09..d48f417 100644 --- a/R/as.Dtable.data.frame.R +++ b/R/as.Dtable.data.frame.R @@ -9,8 +9,8 @@ as.Dtable.data.frame <- function (X, ...) { return ( Dtable( Dmatrix, - PointType = spatstat.geom::marks(X)$PointType, - PointWeight = spatstat.geom::marks(X)$PointWeight + PointType = marks(X)$PointType, + PointWeight = marks(X)$PointWeight ) ) } diff --git a/R/gEnvelope.R b/R/gEnvelope.R index 67dcf2a..e563fc8 100644 --- a/R/gEnvelope.R +++ b/R/gEnvelope.R @@ -1,31 +1,62 @@ -gEnvelope <- -function(X, r = NULL, NumberOfSimulations = 100, Alpha = 0.05, - ReferenceType = "", NeighborType = "", SimulationType = "RandomPosition", - Global = FALSE, verbose = interactive()) { +gEnvelope <- function( + X, + r = NULL, + NumberOfSimulations = 100, + Alpha = 0.05, + ReferenceType = "", + NeighborType = "", + SimulationType = "RandomPosition", + Global = FALSE, + verbose = interactive()) { CheckdbmssArguments() # Choose the null hypothesis - SimulatedPP <- switch (SimulationType, - RandomPosition = expression(rRandomPositionK(X, CheckArguments = FALSE)), - RandomLabeling = expression(rRandomLabeling(X, CheckArguments = FALSE)), - PopulationIndependence = expression(rPopulationIndependenceK(X, ReferenceType, NeighborType, CheckArguments = FALSE)) - ) - if (is.null(SimulatedPP)) - stop(paste("The null hypothesis", sQuote(SimulationType), "has not been recognized.")) + SimulatedPP <- switch ( + SimulationType, + RandomPosition = expression(rRandomPositionK(X, CheckArguments = FALSE)), + RandomLabeling = expression(rRandomLabeling(X, CheckArguments = FALSE)), + PopulationIndependence = expression( + rPopulationIndependenceK( + X, + ReferenceType = ReferenceType, + NeighborType = NeighborType, + CheckArguments = FALSE + ) + ) + ) + if (is.null(SimulatedPP)) { + stop( + paste( + "The null hypothesis", + sQuote(SimulationType), + "has not been recognized." + ) + ) + } + # local envelope, keep extreme values for lo and hi (nrank=1) - Envelope <- envelope(X, fun=ghat, nsim=NumberOfSimulations, nrank=1, - r=r, ReferenceType=ReferenceType, NeighborType=NeighborType, - CheckArguments = FALSE, - simulate=SimulatedPP, verbose=verbose, savefuns=TRUE - ) - attr(Envelope, "einfo")$H0 <- switch (SimulationType, - RandomPosition = "Random Position", - RandomLabeling = "Random Labeling", - PopulationIndependence = "Population Independence" - ) + Envelope <- envelope( + X, + fun = ghat, + nsim = NumberOfSimulations, + nrank = 1, + r = r, + ReferenceType = ReferenceType, + NeighborType = NeighborType, + CheckArguments = FALSE, + simulate = SimulatedPP, + verbose = verbose, + savefuns = TRUE + ) + attr(Envelope, "einfo")$H0 <- switch ( + SimulationType, + RandomPosition = "Random Position", + RandomLabeling = "Random Labeling", + PopulationIndependence = "Population Independence" + ) # Calculate confidence intervals - Envelope <- FillEnvelope(Envelope, Alpha, Global) + Envelope <- FillEnvelope(Envelope, Alpha = Alpha, Global = Global) # Return the envelope return (Envelope) } diff --git a/R/generic.spatstat.R b/R/generic.spatstat.R index 9b26fbd..df69820 100644 --- a/R/generic.spatstat.R +++ b/R/generic.spatstat.R @@ -1,12 +1,10 @@ -sharpen.wmppp <- function(X, ...) -{ +sharpen.wmppp <- function(X, ...) { X <- sharpen.ppp(X, ...) class(X) <- c("wmppp", "ppp") return(X) } -split.wmppp <- function(...) -{ +split.wmppp <- function(...) { Xlist <- split.ppp(...) as.wmppplist <- function(X) { class(X) <- c("wmppp", "ppp") @@ -16,23 +14,20 @@ split.wmppp <- function(...) return(Xlist) } -superimpose.wmppp <- function(...) -{ +superimpose.wmppp <- function(...) { X <- superimpose.ppp(...) class(X) <- c("wmppp", "ppp") return(X) } -unique.wmppp <- function(x, ...) -{ +unique.wmppp <- function(x, ...) { X <- unique.ppp(x, ...) class(X) <- c("wmppp", "ppp") return(X) } -"[.wmppp" <- function(i, j, drop=FALSE, ..., clip=FALSE) -{ - X <- "[.ppp"(i, j, drop=drop, ..., clip=clip) +"[.wmppp" <- function(i, j, drop = FALSE, ..., clip = FALSE) { + X <- "[.ppp"(i, j, drop = drop, ..., clip = clip) class(X) <- c("wmppp", "ppp") return(X) } diff --git a/R/ghat.R b/R/ghat.R index 336e078..5b0f402 100644 --- a/R/ghat.R +++ b/R/ghat.R @@ -1,10 +1,14 @@ -ghat <- -function(X, r = NULL, ReferenceType = "", NeighborType = "", CheckArguments = TRUE) { +ghat <- function( + X, + r = NULL, + ReferenceType = "", + NeighborType = "", + CheckArguments = TRUE) { if (CheckArguments) { CheckdbmssArguments() # Eliminate erroneous configurations - if ((ReferenceType=="" | NeighborType=="") & (ReferenceType!=NeighborType)) { + if ((ReferenceType == "" | NeighborType == "") & (ReferenceType != NeighborType)) { stop("Either two or no point type must be specified.") } } @@ -14,20 +18,20 @@ function(X, r = NULL, ReferenceType = "", NeighborType = "", CheckArguments = TR # Calculate densities # g intra - if (ReferenceType=="" & NeighborType=="") { - lambdaI <- lambdaJ <- X$n/area + if (ReferenceType == "" & NeighborType == "") { + lambdaI <- lambdaJ <- X$n / area } else { # g intra for a single point type - if (ReferenceType==NeighborType) { - X.reduced <- X[marks(X)$PointType==ReferenceType] - lambdaI <- lambdaJ <- X.reduced$n/area + if (ReferenceType == NeighborType) { + X.reduced <- X[marks(X)$PointType == ReferenceType] + lambdaI <- lambdaJ <- X.reduced$n / area } # g inter - if (ReferenceType!=NeighborType) { - X.cross <- X[marks(X)$PointType==ReferenceType] - Y.cross <- X[marks(X)$PointType==NeighborType] - lambdaI <- X.cross$n/area - lambdaJ <- Y.cross$n/area + if (ReferenceType != NeighborType) { + X.cross <- X[marks(X)$PointType == ReferenceType] + Y.cross <- X[marks(X)$PointType == NeighborType] + lambdaI <- X.cross$n / area + lambdaJ <- Y.cross$n / area } } @@ -35,23 +39,25 @@ function(X, r = NULL, ReferenceType = "", NeighborType = "", CheckArguments = TR rmaxdefault <- rmax.rule("K", X$window, lambdaJ) breaks <- handle.r.b.args(window = X$window, rmaxdefault = rmaxdefault) rBest <- breaks$r - denargs <- spatstat.utils::resolve.defaults(list(kernel = "epanechnikov", n = length(rBest), from = 0, to = max(rBest))) + denargs <- spatstat.utils::resolve.defaults( + list(kernel = "epanechnikov", n = length(rBest), from = 0, to = max(rBest)) + ) if (is.null(r)) { r <- rBest } # Find pairs # g intra - if (ReferenceType=="" & NeighborType=="") { - Pairs <- closepairs(X, max(r)) + if (ReferenceType == "" & NeighborType == "") { + Pairs <- closepairs(X, rmax = max(r)) } else { # g intra for a single point type - if (ReferenceType==NeighborType) { - Pairs <- closepairs(X.reduced, max(r)) + if (ReferenceType == NeighborType) { + Pairs <- closepairs(X.reduced, rmax = max(r)) } # g inter if (ReferenceType!=NeighborType) { - Pairs <- crosspairs(X.cross, Y.cross, max(r)) + Pairs <- crosspairs(X = X.cross, Y = Y.cross, rmax = max(r)) } } @@ -63,20 +69,29 @@ function(X, r = NULL, ReferenceType = "", NeighborType = "", CheckArguments = TR # Edge-effect correction if (is.rectangle(X$window) | is.polygonal(X$window)) { - edgewt <- edge.Ripley(XI, matrix(Pairs$d, ncol = 1)) + edgewt <- edge.Ripley(XI, r = matrix(Pairs$d, ncol = 1)) valu <- "iso" desc <- "Ripley isotropic correction estimate of %s" } else { - edgewt <- edge.Trans(XI, XJ, paired = TRUE) + edgewt <- edge.Trans(X = XI, Y = XJ, paired = TRUE) valu <- "trans" desc <- "translation-corrected estimate of %s" } # Estimate g - gEstimate <- sewpcf(Pairs$d, edgewt, denargs, lambdaI*lambdaJ*area) + gEstimate <- sewpcf( + Pairs$d, + w = edgewt, + denargs = denargs, + lambda2area = lambdaI * lambdaJ * area + ) # Calculate values for r if specified if (!autor) { - g <- stats::approx(gEstimate$r, gEstimate$g, xout=r)$y + g <- stats::approx( + x = gEstimate$r, + y = gEstimate$g, + xout = r + )$y gEstimate <- data.frame(r, g) } # Add theoretical value @@ -86,11 +101,18 @@ function(X, r = NULL, ReferenceType = "", NeighborType = "", CheckArguments = TR colnames(gEstimate) <- ColNames # Return the values of g(r) - g <- fv(gEstimate, argu="r", ylab=quote(g(r)), valu=valu, fmla= ". ~ r", - alim=c(0, max(r)), - labl=c("r", "%s[pois](r)", paste("hat(%s)[", valu, "](r)", sep="")), - desc=c("distance argument r", "theoretical Poisson %s", desc), - unitname=X$window$unit, fname="g") + g <- fv( + gEstimate, + argu = "r", + ylab = quote(g(r)), + valu = valu, + fmla= ". ~ r", + alim = c(0, max(r)), + labl = c("r", "%s[pois](r)", paste("hat(%s)[", valu, "](r)", sep="")), + desc = c("distance argument r", "theoretical Poisson %s", desc), + unitname = X$window$unit, + fname = "g" + ) fvnames(g, ".") <- ColNames[-1] return(g) } diff --git a/R/is.wmppp.R b/R/is.wmppp.R index 3afba66..50bb9cf 100644 --- a/R/is.wmppp.R +++ b/R/is.wmppp.R @@ -1,4 +1,3 @@ -is.wmppp <- -function (X) { +is.wmppp <- function (X) { inherits(X, "wmppp") } diff --git a/R/rPopulationIndependenceK.R b/R/rPopulationIndependenceK.R index 920e1ad..96dccb1 100644 --- a/R/rPopulationIndependenceK.R +++ b/R/rPopulationIndependenceK.R @@ -1,20 +1,25 @@ -rPopulationIndependenceK <- -function (X, ReferenceType, NeighborType, CheckArguments = TRUE) { +rPopulationIndependenceK <- function( + X, + ReferenceType, + NeighborType, + CheckArguments = TRUE) { - if (CheckArguments) + if (CheckArguments) { CheckdbmssArguments() - + } + # Eliminate useless points - X.reduced <- X[marks(X)$PointType==ReferenceType | marks(X)$PointType==NeighborType] + X.reduced <- + X[marks(X)$PointType == ReferenceType | marks(X)$PointType == NeighborType] RandomizedX <- X.reduced # Reduce the factor levels to two (factor eliminates the levels with no points) - Marks <- factor(spatstat.geom::marks(X.reduced)$PointType) + Marks <- factor(marks(X.reduced)$PointType) # The new point pattern has classical spatstat marks RandomizedX <- RandomizedX %mark% Marks # Split reference and neighbor points X.split <- split(RandomizedX) # Randomly shift the neighbors - rshift(X.split, which=NeighborType) -> RandomizedX.split + rshift(X.split, which = NeighborType) -> RandomizedX.split # Reunify the split point pattern RandomizedX.split -> split(RandomizedX) # Reorganize the marks (add weight) diff --git a/R/rRandomPositionK.R b/R/rRandomPositionK.R index 73a5b7c..2c0c74e 100644 --- a/R/rRandomPositionK.R +++ b/R/rRandomPositionK.R @@ -1,11 +1,14 @@ -rRandomPositionK <- -function(X, Precision = 0, CheckArguments = TRUE) { +rRandomPositionK <- function( + X, + Precision = 0, + CheckArguments = TRUE) { - if (CheckArguments) + if (CheckArguments) { CheckdbmssArguments() - + } + # Draw in a binomial process - RandomizedX <- runifpoint(X$n, win=X$window) + RandomizedX <- runifpoint(X$n, win = X$window) # Precision if (Precision) { @@ -14,7 +17,10 @@ function(X, Precision = 0, CheckArguments = TRUE) { } # Apply original marks to new points - marks(RandomizedX) <- data.frame(PointWeight=marks(X)$PointWeight, PointType=spatstat.geom::marks(X)$PointType) + marks(RandomizedX) <- data.frame( + PointWeight = marks(X)$PointWeight, + PointType = marks(X)$PointType + ) class(RandomizedX) <- c("wmppp", "ppp") return (RandomizedX) }