sdr.R
#'
#' sdr.R
#'
#' Sufficient Dimension Reduction
#'
#' Matlab original: Yongtao Guan
#' Translated to R by: Suman Rakshit
#' Adapted for spatstat: Adrian Baddeley
#'
#' GNU Public Licence 2.0 || 3.0
#'
#' $Revision: 1.15 $ $Date: 2020/01/30 05:10:49 $
#'
sdr <- function(X, covariates, ...) {
UseMethod("sdr")
}
sdr.ppp <- local({
sdr.ppp <- function(X, covariates,
method=c("DR", "NNIR", "SAVE", "SIR", "TSE"),
Dim1=1, Dim2=1, predict=FALSE, ...) {
stopifnot(is.ppp(X))
method <- match.arg(method)
trap.extra.arguments(...)
#' ensure 'covariates' is a list of compatible images
if(!inherits(covariates, "imlist") && !all(sapply(covariates, is.im)))
stop("Argument 'covariates' must be a list of images")
nc <- length(covariates)
if(nc == 0)
stop("Need at least one covariate!")
if(nc < Dim1 + (method == "TSE") * Dim2)
stop(paste(if(method == "TSE") "Dim1 + Dim2" else "Dim1",
"must not exceed the number of covariates"),
call.=FALSE)
if(nc > 1 && !do.call(compatible, unname(covariates)))
covariates <- do.call(harmonise, covariates)
#' extract corresponding pixel values including NA's
Ypixval <- sapply(lapply(covariates, as.matrix), as.vector)
#' compute sample mean and covariance matrix
m <- colMeans(Ypixval, na.rm=TRUE)
V <- cov(Ypixval, use="complete")
#' evaluate each image at point data locations
YX <- sapply(covariates, safelook, Y=X)
#' apply precomputed standardisation
Zx <- t(t(YX) - m) %*% matrixinvsqrt(V)
#' ready
coordsX <- coords(X)
result <-
switch(method,
DR = calc.DR(COV=V, z=Zx, Dim=Dim1),
NNIR = calc.NNIR(COV=V, z=Zx, pos=coordsX, Dim=Dim1),
SAVE = calc.SAVE(COV=V, z=Zx, Dim=Dim1),
SIR = calc.SIR(COV=V, z=Zx ),
TSE = calc.TSE(COV=V, z=Zx, pos=coordsX, Dim1=Dim1, Dim2=Dim2)
)
#'
covnames <- names(covariates) %orifnull% paste0("Y", 1:nc)
dimnames(result$B) <- list(covnames, paste0("B", 1:ncol(result$B)))
if(method == "TSE") {
result$M1 <- namez(result$M1)
result$M2 <- namez(result$M2)
} else {
result$M <- namez(result$M)
}
if(predict) result$Y <- sdrPredict(covariates, result$B)
return(result)
}
safelook <- function(Z, Y, ...) { safelookup(Z, Y, ...) }
namez <- function(M, prefix="Z") {
dimnames(M) <- list(paste0(prefix, 1:nrow(M)),
paste0(prefix, 1:ncol(M)))
return(M)
}
sdr.ppp
})
sdrPredict <- function(covariates, B) {
if(!is.matrix(B)) {
if(is.list(B) && is.matrix(BB <- B$B)) B <- BB else
stop("B should be a matrix, or the result of a call to sdr()",
call.=FALSE)
}
if(!inherits(covariates, "imlist") && !all(sapply(covariates, is.im)))
stop("Argument 'covariates' must be a list of images")
stopifnot(nrow(B) == length(covariates))
result <- vector(mode="list", length=ncol(B))
for(j in seq_along(result)) {
cj <- as.list(B[,j])
Zj <- mapply("*", cj, covariates, SIMPLIFY=FALSE)
result[[j]] <- im.apply(Zj, sum)
}
names(result) <- colnames(B)
return(as.solist(result))
}
##............ DR (Directional Regression) ..........................
calc.DR <- function(COV, z, Dim){
## Description: Naive Directional Regression Method
## Input:
## COV - cov{X(s)}
## z - standardized X(s) on SPP locations
## Dim - the CS dimension
## Output:
## B - the estimated CS basis
## M - the kernel matrix
ss <- nrow(z)
ncov <- ncol(z)
## M1 <- (t(z) %*% z)/ss - diag(1,ncov)
M1 <- crossprod(z)/ss - diag(1,ncov)
M1 <- M1 %*% M1 # the SAVE kernel
covMean <- matrix(colMeans(z),ncol=1)
M2 <- covMean %*% t(covMean)
M3 <- M2 * (base::norm(covMean, type="2"))^2 # the SIR kernel
M2 <- M2 %*% M2 # the SIR-2 kernel
M <- (M1 + M2 + M3)/3 # the DR kernel
SVD <- svd(M)
B <- SVD$u[,1:Dim]
B <- matrixinvsqrt(COV) %*% B # back to original scale
return(list(B=B, M=M))
}
## ............ NNIR (Nearest Neighbor Inverse Regression) ...........
calc.NNIR <- function(COV, z, pos, Dim) {
## Description: Nearest Neighbor Inverse Regression
## Input:
## COV - cov{X(s)}
## z - standardized X(s) on SPP locations
## pos - the position of SPP events
## Dim - the CS dimension
## Output:
## B - the estimated CS basis
## M - the kernel matrix
ss <- nrow(z) # sample size
# ncov <- ncol(z) # predictor dimension
jj <- nnwhich(pos) # identify nearest neighbour of each point
dir <- z - z[jj, , drop=FALSE] # empirical direction
IM <- sumouter(dir) # inverse of kernel matrix: sum of outer(dir[i,], dir[i,])
M <- solve(IM/ss) # invert kernel matrix
SVD <- svd(M)
B <- matrixinvsqrt(COV) %*% SVD$u[, 1:Dim, drop=FALSE]
return(list(B=B, M=M))
}
## ........... SAVE (Sliced Average Variance Estimation) ...........
calc.SAVE <- function(COV, z, Dim){
## Description: Naive Directional Regression Method
## Input
## COV - cov{X(s)}
## z - standardized X(s) on SPP locations
## Dim - the central space dimension
## Value
## B - the estimated CS basis
## M - the kernel matrix
# ss <- nrow(z)
ncov <- ncol(z)
M <- diag(1,ncov) - cov(z)
M <- M %*% M
SVD <- svd(M)
B <- SVD$u[,1:Dim]
B <- matrixinvsqrt(COV) %*% B
return(list(B=B, M=M))
}
##.......... SIR (Sliced Inverse Regression) ......................
calc.SIR <- function(COV, z){
## Description: Naive Directional Regression Method
## Input:
## COV - cov{X(s)}
## z - standardized X(s) on SPP locations
## Output:
## B - the estimated CS basis
## M - the kernel matrix
covMean <- colMeans(z)
B <- matrixinvsqrt(COV) %*% covMean # do SIR estimation
B <- B/sqrt(sum(B^2)) # normalise to unit length
M <- covMean %*% t(covMean) # create kernel matrix
return(list(B=B, M=M))
}
## ............. TSE (Two-Step Estimation) ....................
calc.TSE <- function(COV, z, pos, Dim1, Dim2) {
## Description: A Two-Step Method
## Input:
## COV - cov{X(s)}
## z - standardized X(s) on SPP locations
## Dim1 - the S1 dimension
## Dim2 - the S2 dimension
## Output:
## B - the estimated CS basis. Its first Dim1 columns
## are estimating S1 and the remaining Dim2 columns are
## estimating S2. In case of null space, a zero vector is reported.
## M1 - the kernel matrix of DR
## M2 - the kernel matrix of NNIR, which might be subject
## to some change, depending on the results of M1.
# ss <- nrow(z) # sample size
ncov <- ncol(z) # predictor dimension
est1 <- calc.DR(COV, z, ncov) # do DR estimation
est2 <- calc.NNIR(COV, z, pos, ncov) # do NNIR estimation
M1 <- est1$M
M2 <- est2$M
if(Dim1 > 0) {
U <- svd(M1)$u
B1 <- U[ , 1:Dim1, drop=FALSE] # get S1 estimate
## Q <- diag(1, ncov) - B1 %*% solve(t(B1) %*% B1) %*% t(B1)
Q <- diag(1, ncov) - B1 %*% solve(crossprod(B1)) %*% t(B1)
# contract orthogonal basis
M2 <- Q %*% M2 %*% Q # do constrained NNIR
} else {
B1 <- matrix(0, ncov, 1)
}
if(Dim2 > 0) {
U <- svd(M2)$u # do SVD for possibly updated M2
B2 <- U[ , 1:Dim2, drop=FALSE] # get basis estimator
} else {
B2 <- matrix(0, ncov, 1)
}
B <- matrixinvsqrt(COV) %*% cbind(B1,B2)
return(list(B=B, M1=M1, M2=M2))
}
## ////////////////// ADDITIONAL FUNCTIONS /////////////////////
subspaceDistance <- function(B0,B1) {
## ======================================================== #
## Evaluate the distance between the two linear spaces S(B0) and S(B1).
## The measure used is the one proposed by Li et al. (2004).
## ======================================================== #
stopifnot(is.matrix(B0))
stopifnot(is.matrix(B1))
## Proj0 <- B0 %*% solve((t(B0) %*% B0)) %*% t(B0) # Proj matrix on S(B0)
Proj0 <- B0 %*% solve(crossprod(B0)) %*% t(B0) # Proj matrix on S(B0)
lam <- svd(B1) # check whether B1 is singular
U <- lam$u
D <- lam$d
# V <- lam$v
B2 <- U[, D > 1e-09] # keep non-singular directions
Proj1 <- B2 %*% solve((t(B2) %*% B2)) %*% t(B2) # Proj matrix on S(B.hat)
Svd <- svd(Proj0 - Proj1) # Do svd for P0-P1
dist <- max(abs(Svd$d)) # Get the maximum absolute svd value
return(dist)
}
dimhat <- function(M){
#' Description: Maximum Descent Estimator for CS Dim
#' Input:
#' M - the estimated kernel matrix
#' Output:
#' dimhat - the estimated CS dim (assume dim>0)
stopifnot(is.matrix(M))
ncov <- ncol(M) # predictor dimension
maxdim <- max((ncov-1), 5) # maximum structure dimension
SVD <- svd(M) # svd of kernel matrix
lam <- SVD$d
eps <- 1e-06
lam <- lam + rep(eps,ncov) # add ridge effect
lam1 <- lam[-ncov]
lam2 <- lam[-1]
dif <- lam1/lam2
dif <- dif[1 : maxdim] # the magnitude of drop
retval <- which.max(dif) # find Maximum Descent estimator
return(retval)
}