R/scgwr_p.R

In scgwr: Scalable Geographically Weighted Regression

scgwr_p    <- function (coords, y, x=NULL, knn = 100, kernel = "gau", p = 4,
                        approach="CV", nsamp=NULL, cl=NULL){

  ic_score <- function(par,X,y,p,MM,MM2,mm,nx,id_x1,id_x2,id_xy,id_p,id_py,nsamp,XX0,Xy0) {
    par2 <- par^2
    R0   <- rep(par2[1], p + 1)
    for (r0 in 2:(p + 1)) R0[r0] <- par2[1]^r0
    R0   <-R0/sum( R0 )

    MMM  <- R0[id_p] * MM
    mmm  <- R0[id_py]* mm
    MMM2 <- 2*par2[2]*( R0[id_p] * MM ) + ( R0[id_p]^2 * MM2 )

    MMMsum <- NULL
    MMMsum2<- NULL
    mmmsum <- NULL
    for (ll1 in 1:nx) {
      for (ll2 in 1:nx) {
        MMM_sub <- colSums(MMM[(id_x1 == ll1) & (id_x2 == ll2), ])
        MMMsum <- rbind(MMMsum, MMM_sub)

        MMM_sub2<- colSums(MMM2[(id_x1 == ll1) & (id_x2 == ll2), ])
        MMMsum2 <- rbind(MMMsum2, MMM_sub2)
      }
      mmm_sub <- colSums(mmm[id_xy == ll1, ])
      mmmsum <- rbind(mmmsum, mmm_sub)
    }
    MMMlist <- list(NULL)
    MMMlist2<- list(NULL)
    mmmlist <- list(NULL)
    xxxlist <- list(NULL)# code is inefficient (it must not needed)
    for (ll3 in 1:nsamp) {
      MMMlist[[ll3]] <- matrix(MMMsum[, ll3], nx, nx) + par2[2] * XX0
      MMMlist2[[ll3]]<- matrix(MMMsum2[, ll3], nx, nx) + par2[2] ^ 2 * XX0
      mmmlist[[ll3]] <- mmmsum[, ll3] + par2[2] * Xy0
      xxxlist[[ll3]] <- X[ll3,]
    }
    tryres <- try(MMMinv <- lapply(MMMlist, solve), silent = TRUE)
    if (class(tryres) == "try-error") {
      obj <- 10^10
    } else {
      beta <- t(matrix(unlist(mapply("%*%", MMMinv, mmmlist,
                                     SIMPLIFY = FALSE)), nrow = nx, ncol = nsamp))

      trS00  <- mapply("%*%", MMMinv , xxxlist,SIMPLIFY = FALSE)
      trS0   <- mapply("%*%", xxxlist, trS00   ,SIMPLIFY = FALSE)
      trS    <- Reduce("+", trS0)#*g00

      trSS00  <- mapply("%*%", MMMlist2 , trS00,   SIMPLIFY = FALSE)
      trSS0   <- mapply(function(x,y) sum(x*y), trS00, trSS00, SIMPLIFY = FALSE)
      trSS    <- Reduce("+", trSS0)

      pred    <- rowSums(X * beta)#[nlist,]
      sse     <- sum((y - pred)^2)#[nlist]
      sig     <- sqrt( sse/nsamp )
      AICc    <- 2*nsamp*log(sig)+nsamp*log(2*pi)+nsamp*(nsamp+trS)/(nsamp-2-trS)
    }

    return(AICc)
  }

  cv_score <- function(par,X,y,p,MM,mm,nx,id_x1,id_x2,id_xy,id_p,id_py,nsamp,XX0,Xy0){
    par2 <- par^2
    R0   <- rep(par2[1], p + 1)
    for (r0 in 2:(p + 1)) R0[r0] <- par2[1]^r0
    R0   <-R0/sum( R0 )

    MMM  <- R0[id_p] * MM#[,nlist]
    mmm  <- R0[id_py] * mm#[,nlist]
    MMMsum <- NULL
    mmmsum <- NULL
    for (ll1 in 1:nx) {
      for (ll2 in 1:nx) {
        MMM_sub<- colSums(MMM[(id_x1 == ll1) & (id_x2 == ll2), ])
        MMMsum <- rbind(MMMsum, MMM_sub)
      }
      mmm_sub <- colSums(mmm[id_xy == ll1, ])
      mmmsum <- rbind(mmmsum, mmm_sub)
    }
    MMMlist <- list(NULL)
    mmmlist <- list(NULL)

    for (ll3 in 1:nsamp) {
      MMMlist[[ll3]] <- matrix(MMMsum[, ll3], nx, nx) + par2[2] * XX0
      mmmlist[[ll3]] <- mmmsum[, ll3] + par2[2] * Xy0
    }
    tryres <- try(MMMinv <- lapply(MMMlist, solve), silent = TRUE)
    if (class(tryres) == "try-error") {
      error <- 10^10
    }
    else {
      beta <- t(matrix(unlist(mapply("%*%", MMMinv, mmmlist,
                                     SIMPLIFY = FALSE)), nrow = nx, ncol = nsamp))
      pred <- rowSums(as.matrix(X) * beta)#[nlist,]
      error <- sum((y - pred)^2)#[nlist]
    }
    return(error)
  }

  proc  <-function(coords, X, y, coords_samp, X_samp, y_samp, knn, approach, kernel,p){
    n_samp <- length(y_samp)
    Dnn    <- get.knnx(coords,coords_samp, knn-1)
    if( approach != "CV"){
      Dnn$nn.index<-as.matrix( cbind( 1:n_samp, Dnn$nn.index) )
      Dnn$nn.dist <-as.matrix( cbind( 0  , Dnn$nn.dist ) )
    }

    bmax <- sqrt((max(coords[,1])-min(coords[,1]))^2+(max(coords[,2])-min(coords[,2]))^2)
    Dband<- seq(0,bmax,len=10000)
    Did  <- Dnn$nn.index
    if (kernel == "gau") {
      ban0  <- median(Dnn$nn.dist[, min(50,knn)])/sqrt(3)
      G0    <- exp(-(Dnn$nn.dist/ban0)^2)
      G0band<- exp(-(Dband/ban0)^2)
    } else if (kernel == "exp") {
      ban0  <- median(Dnn$nn.dist[, min(50,knn)])/3
      G0    <- exp(-Dnn$nn.dist/ban0)
      G0band<- exp(-Dband/ban0)
    }

    XX0   <- crossprod(X)
    Xy0   <- crossprod(X,y)
    id_p  <- rep(rep(1:(p + 1), each = nx), nx)
    id_x1 <- rep(1:nx, each = nx * (p + 1))
    id_x2 <- rep(1:nx, nx * (p + 1))
    id_py <- id_p [id_x1 == 1]
    id_xy <- id_x2[id_x1 == 1]
    MM    <- matrix(0, ncol = n_samp, nrow = length(id_p))
    mm    <- matrix(0, ncol = n_samp, nrow = sum(id_x1 == 1))
    if(approach=="AICc"){
      MM2 <- matrix(0, ncol = n_samp, nrow = length(id_p))
      G   <- matrix(1, nrow = knn, ncol = p + 1)
    } else {
      MM2 <- NULL
      G   <- matrix(1, nrow = knn - 1, ncol = p + 1)
    }

    for (i in 1:n_samp){
      for (p0 in 1:p) G[, p0 + 1] <- G0[i, ]^(2^(p/2)/2^p0)
      X_sub <- as.matrix(X[Did[i, ],])
      y_sub <- y[Did[i, ]]
      kk0   <- 1
      for (k0 in 1:nx) {
        GXp   <- G * X_sub[, k0]
        if(approach=="AICc"){
          GXp2<- G^2 * X_sub[, k0]
        }
        for (k1 in 1:(p + 1)) {
          MM [(id_x1 == k0) & (id_p == k1), i]  <- colSums(GXp [,k1] * X_sub)
          mm [(id_xy == k0) & (id_py == k1), i] <- sum(GXp[,k1] * y_sub)
          if(approach=="AICc"){
            MM2[(id_x1 == k0) & (id_p == k1), i]  <- colSums(GXp2[,k1] * X_sub)
          }
        }
      }
    }

    return(list(MM=MM, MM2=MM2, mm=mm, XX0=XX0,Xy0=Xy0,id_p=id_p,id_x1=id_x1,id_x2=id_x2,
                id_py=id_py,id_xy=id_xy,G0=G0,Dnn=Dnn, ban0=ban0))
  }

  ##########################################################


  n <- length(y)
  if (knn >= n) knn   <- n - 1
  if (is.null(x)) {
    X     <- as.matrix(rep(1, n))
    xname <- "(Intercept)"
    x_id <- NULL
  } else {
    X00 <- as.matrix(x)
    if (is.numeric(X00) == F) {
      mode(X00) <- "numeric"
    }

    if( dim( as.matrix(x) )[2]==1 ){
      x_id <- (sd(x) !=0)
    } else {
      x_id <- (apply(x,2,sd) !=0)#(x %>% summarise_if(is.numeric, funs(sd)) !=0)
    }
    x_nm <- colnames(x)[x_id] #added

    if( all(x_id == 0)){
      X <- as.matrix(rep(1, n))
      xname <- "(Intercept)"
      x_id <- NULL
    } else {
      X <- as.matrix(cbind(1, X00[, x_id]))
      xname <- c("(Intercept)", names(as.data.frame( X00[, x_id])))
    }
  }
  nx  <- dim(X)[2]

  if(approach=="AICc"){
    if( !is.null( nsamp ) ){
      message("NOTE: nsamp is ignored because it is not available for the AICc-based estimation")
    }

    nlist<-1:n
    nsamp<-n
  } else if( !is.null(nsamp) ){
    if(nsamp < n){
      nlist<-sample(n,nsamp)
    } else {
      nlist<-1:n
      nsamp<-n
    }
  } else {
    nlist<-1:n
    nsamp<-n
  }

  coords_samp<- coords[nlist,]
  X_samp     <- X[nlist,]
  y_samp     <- y[nlist]
  misc       <- proc(coords=coords,coords_samp=coords_samp, X=X, y=y, X_samp=X_samp, y_samp=y_samp, knn=knn, approach=approach, kernel=kernel, p=p)
  MM         <- misc$MM
  MM2        <- misc$MM2
  mm         <- misc$mm
  XX0        <- misc$XX0
  Xy0        <- misc$Xy0
  G0         <- misc$G0
  ban0       <- misc$ban0
  id_p       <- misc$id_p
  id_x1      <- misc$id_x1
  id_x2      <- misc$id_x2
  id_py      <- misc$id_py
  id_xy      <- misc$id_xy


  if(is.null(cl)) {
    cl <- makeCluster(detectCores(),setup_strategy = "sequential")
  } else {
    cl <- makeCluster(cl,setup_strategy = "sequential")
  }
  setDefaultCluster(cl=cl)

  par0       <- c(1, 0.01)
  if(approach == "CV"){
    res0    <- optimParallel(par=par0,fn = cv_score, X=X_samp,y=y_samp,p=p, MM=MM,mm=mm,nx=nx,id_x1=id_x1,
                             id_x2=id_x2,id_xy=id_xy,id_p=id_p,id_py=id_py,nsamp=nsamp,XX0=XX0,Xy0=Xy0,
                             lower=c(-10,-Inf),upper=c(10,Inf),parallel=list(forward=FALSE,loginfo = FALSE))
    cvscore <- sqrt(res0$value/n)
  } else {
    res0    <- optimParallel(par=par0,fn = ic_score, X=X_samp,y=y_samp,p=p,MM=MM,MM2=MM2,mm=mm,nx=nx,id_x1=id_x1,
                             id_x2=id_x2,id_xy=id_xy,id_p=id_p,id_py=id_py,nsamp=nsamp,XX0=XX0,Xy0=Xy0,
                             lower=c(-10,-Inf),upper=c(10,Inf),parallel=list(forward=FALSE,loginfo = FALSE))
    cvscore <-NULL
  }
  setDefaultCluster(cl=NULL); stopCluster(cl)

  par    <- res0$par
  par2   <- par^2
  if( 100001 >= n){#nsamp == n
    niter<-1
    ilist<-list( 1:n )
  } else {
    niter<-round(n/100000)
    ilist<-list( 1:100000 )
    nmax <-0
    ii   <-1
    while(nmax<=n){
      ilist0     <- c(nmax + 1:100000)
      ilist0_id  <- ilist0 <= n
      if( sum( ilist0_id ) >= 1 ){
        ilist[[ii]]<- ilist0[ilist0_id]
      }
      nmax<-nmax+100000
      ii  <-ii+1
    }
  }

  BETA       <-matrix(0,nrow=n,ncol=nx)
  BSE_no_sig <-BETA
  ENP        <-0
  PRED       <-rep(0,n)
  for(iii in 1:length(ilist)){
    jd  <- ilist[[iii]]
    n_iii<-length(jd)
    Dnn <- get.knnx(coords,coords[jd,], knn)#not knn -1 because of using get.knnx
    Did  <- Dnn$nn.index
    if (kernel == "gau") {
      G0    <- exp(-(Dnn$nn.dist/ban0)^2)
    } else if (kernel == "exp") {
      G0    <- exp(-Dnn$nn.dist/ban0)
    }

    MM <- matrix(0, ncol = n_iii, nrow = length(id_p))
    MM2<- matrix(0, ncol = n_iii, nrow = length(id_p))
    mm <- matrix(0, ncol = n_iii, nrow = sum(id_x1 == 1))
    for (i in 1:n_iii) {
      G <- matrix(1, nrow = knn, ncol = p + 1)
      for (p0 in 1:p) G[, p0 + 1] <- G0[i, ]^(2^(p/2)/2^p0)
      X_sub <- as.matrix(X[Did[i, ],])
      y_sub <- y[Did[i, ]]
      kk0   <- 1
      for (k0 in 1:nx) {
        GXp <- G   * X_sub[, k0]
        GXp2<- G^2 * X_sub[, k0]
        for (k1 in 1:(p + 1)) {
          MM [(id_x1 == k0) & (id_p == k1), i]  <- colSums(GXp [,k1] * X_sub)
          MM2[(id_x1 == k0) & (id_p == k1), i]  <- colSums(GXp2[,k1] * X_sub)
          mm [(id_xy == k0) & (id_py == k1), i] <- sum(GXp[,k1] * y_sub)
        }
      }
    }

    R0   <- rep(par2[1], p + 1)
    for (r0 in 2:(p + 1)) R0[r0] <- par2[1]^r0
    R0   <-R0/sum( R0 )

    MMM  <- R0[id_p] * MM
    mmm  <- R0[id_py] * mm
    MMM2 <- 2*par2[2]*MMM + R0[id_p]^2 * MM2

    MMMsum <- NULL
    MMMsum2<- NULL
    mmmsum <- NULL
    for (ll1 in 1:nx) {
      for (ll2 in 1:nx) {
        MMM_sub <- colSums(MMM[(id_x1 == ll1) & (id_x2 == ll2), ])
        MMMsum  <- rbind(MMMsum, MMM_sub)

        MMM_sub2<- colSums(MMM2[(id_x1 == ll1) & (id_x2 == ll2), ])
        MMMsum2 <- rbind(MMMsum2, MMM_sub2)
      }
      mmm_sub <- colSums(mmm[id_xy == ll1, ])
      mmmsum <- rbind(mmmsum, mmm_sub)
    }

    MMMlist <- list(NULL)
    MMMlist2<- list(NULL)
    mmmlist <- list(NULL)
    xxxlist <- list(NULL)# code is inefficient (it must not needed)
    for (ll3 in 1:n_iii) {
      MMMlist[[ll3]] <- matrix(MMMsum[, ll3], nx, nx)  + par2[2] * XX0
      MMMlist2[[ll3]]<- matrix(MMMsum2[, ll3], nx, nx) + par2[2] ^ 2 * XX0
      mmmlist[[ll3]] <- mmmsum[, ll3] + par2[2] * Xy0
      xxxlist[[ll3]] <- as.matrix( X[ilist[[iii]],] )[ll3,]
    }

    tryres  <- try(MMMinv <- lapply(MMMlist, solve), silent = TRUE)
    beta    <- t(matrix(unlist(mapply("%*%", MMMinv, mmmlist,
                                      SIMPLIFY = FALSE)), nrow = nx, ncol = n_iii))

    trS00  <- mapply("%*%", MMMinv , xxxlist,SIMPLIFY = FALSE)
    trS0   <- mapply("%*%", xxxlist, trS00   ,SIMPLIFY = FALSE)
    trS    <- Reduce("+", trS0)#*g00

    trSS00  <- mapply("%*%", MMMlist2 , trS00,   SIMPLIFY = FALSE)
    trSS0   <- mapply(function(x,y) sum(x*y), trS00, trSS00, SIMPLIFY = FALSE)
    trSS    <- Reduce("+", trSS0)
    enp     <- 2*trS - trSS

    pred    <- rowSums(X[ilist[[iii]],] * beta)
    bse00   <- mapply("%*%", MMMlist2, MMMinv,SIMPLIFY = FALSE)
    bse0    <- mapply("%*%", MMMinv  , bse00 ,SIMPLIFY = FALSE)
    bse0    <- lapply(bse0, function(x) sqrt(diag( x )))
    bse     <- matrix(unlist(bse0), nrow = n_iii, ncol = nx, byrow=TRUE)

    BETA[ilist[[iii]],]       <- beta
    BSE_no_sig[ilist[[iii]],] <- bse
    ENP     <-ENP + enp
    PRED[ilist[[iii]]]        <- pred
  }

  pred    <- PRED
  resid   <- y - pred
  sse     <- sum( resid ^ 2 )
  sig     <- sqrt( sse/( n - enp ) )
  bse     <- c(sig)*BSE_no_sig#sqrt(sig) is replaced with sig
  R2      <- cor(pred, y)^2
  adjR2   <- 1 - (1 - R2) * (n - 1)/(n - enp - 1)

  beta    <- BETA
  bt      <- beta/bse
  bp0     <- data.frame(2 - 2 * pt(abs(bt), df = c(n) - enp))
  bp      <- bp0    #bp<- data.frame(bp0 * c(abs(enp/nx))) # this part is changed
  bp[bp > 1] <- 1

  beta    <- data.frame(beta)
  bse     <- data.frame(bse)
  bt      <- data.frame(beta/bse)
  names(beta) <- names(bse) <- names(bt) <- names(bp) <-names(bp0) <- xname

  r       <- par2[1]
  pen     <- par2[2]
  spar    <- data.frame(c(r, pen))
  names(spar)    <- "Estimate"
  rownames(spar) <- c("scale", "penalty")

  sig0    <- sqrt( sse/n )
  loglik  <-  -n*log(sig0)- n/2*log(2*pi)
  AICc    <- -2*loglik + n*(n+trS)/(n-2-trS)

  if(approach=="CV"){
    e_stat <- data.frame(stat = c( sig, R2, adjR2, loglik,AICc,cvscore))#sqrt(sig) is replaced with sig
    rownames(e_stat) <- c("resid_SE", "R2", "adjR2", "logLik","AICc","CV_score(RMSE)")
  } else {
    e_stat <- data.frame(stat = c( sig, R2, adjR2, loglik,AICc))   #sqrt(sig) is replaced with sig
    rownames(e_stat) <- c("resid_SE", "R2", "adjR2", "logLik","AICc")
  }

  other <- list( kernel = kernel, knn = knn, x_id = x_id, xname = xname, ban0 = ban0, nx=nx,
                 p=p, X=X, y=y, par = par, XX0 = XX0, Xy0 = Xy0, enp = enp, coords=coords )

  result<- list(b = beta, bse = bse, t = bt, p = bp0, par = spar,
                e = e_stat, pred = pred, resid = resid, other = other, call = match.call() )
  class( result ) <- "scgwr_p"
  return( result )
}

print.scgwr_p <- function(x, ...)
{
  cat("Call:\n")
  print(x$call)
  cat("\n----Spatially varying coefficients on x (summary)----\n")
  cat("\nCoefficient estimates:\n")
  print( summary( x$b ) )
  cat("\nStatistical significance:\n")
  p01<-apply(x$p,2,function(x) sum(x<0.01))
  p05<-apply(x$p,2,function(x) sum(x<0.05)) - p01
  p10<-apply(x$p,2,function(x) sum(x<0.10)) - p01 - p05
  p90<-length(x$p[,1]) - p01 - p05 - p10
  pv <-data.frame(rbind( p90, p10, p05, p01))
  names(pv)[1]  <- "Intercept"
  row.names(pv) <- c("Not significant", "Significant (10% level)",
                     "Significant ( 5% level)","Significant ( 1% level)")
  print(pv)
  cat("\n----Variance parameters----------------------------------\n")
  print( x$par )
  cat("\n----Error statistics-------------------------------------\n")
  print(x$e)
  invisible(x)
}