R/MGWRSAR.R

Defines functions MGWRSAR

Documented in MGWRSAR

#' Estimation of linear and local linear model with spatial autocorrelation model (mgwrsar).
#'
#' MGWRSAR is is a wrapper function for estimating linear and local linear models
#' with spatial autocorrelation (SAR models with spatially varying coefficients).
#'
#' @usage MGWRSAR(formula,data,coord,fixed_vars=NULL,kernels,H,
#' Model='GWR',control=list())
#'
#' @param formula  a formula.
#' @param fixed_vars a vector with the names of spatiallay constant coefficient for
#' mixed model. All other variables present in formula are supposed to be spatially
#' varying. If empty or NULL (default), all variables in formula are supposed to be
#' spatially varying.
#' @param data a dataframe or a spatial dataframe (sp package).
#' @param coord default NULL, a dataframe or a matrix with coordinates, not
#' required if data is a spatial dataframe.
#' @param Model character containing the type of model: Possible values are "OLS",
#' "SAR", "GWR" (default), "MGWR" , "MGWRSAR_0_0_kv","MGWRSAR_1_0_kv",
#' "MGWRSAR_0_kc_kv", "MGWRSAR_1_kc_kv", "MGWRSAR_1_kc_0". See Details for more
#' explanation.
#' @param kernels A vector containing the kernel types. Possible types:
#' rectangle ("rectangle"), bisquare ("bisq"), tricube ("tcub"), epanechnikov ("epane"), gaussian
#' ("gauss")) .
#' @param H vector containing the bandwidth parameters for the kernel functions.
#' @param control list of extra control arguments for MGWRSAR wrapper - see Details below
#' @details
#' \itemize{
#' \item{Z}{ a matrix of variables for genralized kernel product, default NULL.}
#' \item{W}{ a row-standardized spatial weight matrix for Spatial Aurocorrelation, default NULL.}
#' \item{type}{ verbose mode, default FALSE.}
#' \item{adaptive}{A vector of boolean to choose adaptive version for each kernel.}
#' \item{kernel_w}{ the type of kernel for computing W, default NULL.}
#' \item{h_w}{ the bandwidth value for computing W, default 0.}
#' \item{Method}{ estimation technique for computing the models with Spatial Dependence. '2SLS' or 'B2SLS', default '2SLS'.}
#' \item{TP}{ A vector of target points, default NULL.}
#' \item{isgcv}{ computing LOOCV criteria (for example for selecting optimal bandwidth), default FALSE.}
#' \item{isfgcv}{ if TRUE, simplify the computation of CV criteria (remove or not i
#' when using local instruments for model with lambda spatially varying), default TRUE.}
#' \item{maxknn}{ when n >NmaxDist, only the maxknn first neighbours are used for distance compution, default 500.}
#' \item{NmaxDist}{ when n >NmaxDist only the maxknn first neighbours are used for distance compution, default 5000}
#' \item{verbose}{ verbose mode, default FALSE.}
#'}
#' @return MGWRSAR returns an object  of class mgwrsar with at least the following components:
#' \item{Betav}{ matrix of coefficients of dim(n,kv) x kv.}
#' \item{Betac}{ vector of coefficients of length kc.}
#' \item{Model}{ The sum of square residuals.}
#' \item{Y}{ The dependent variable.}
#' \item{XC}{ The explanatory variables with constant coefficients.}
#' \item{XV}{ The explanatory variables with varying coefficients.}
#' \item{X}{ The explanatory variables.}
#' \item{W}{ The spatial weight matrix for spatial dependence.}
#' \item{isgcv}{ if gcv has been computed.}
#' \item{edf}{ The estimated degrees of freedom.}
#' \item{formula}{The formula.}
#' \item{data}{ The dataframe used for computation.}
#' \item{Method}{ The type of model.}
#' \item{coord}{ The spatial coordinates of observations.}
#' \item{H}{ The bandwidth vector.}
#' \item{fixed_vars}{ The names of constant coefficients.}
#' \item{kernels}{ The kernel vector.}
#' \item{SSR}{ The sum of square residuals.}
#' \item{residuals}{ The vector of residuals.}
#' \item{fit}{ the vector of fitted values.}
#' \item{sev}{ local standard error of parameters.}
#' \item{NN}{ Maximum number of neighbors for weights computation}
#'
#' MGWRSAR is is a wrapper function for estimating linear and local linear model
#' with spatial autocorrelation that  allows to estimate the following models :
#' \eqn{y=\beta_c X_c+\,\epsilon_i} (OLS)
#'
#' \eqn{y=\beta_v(u_i,v_i) X_v+\,\epsilon_i} (GWR)
#'
#' \eqn{y=\beta_c X_c+\beta_v(u_i,v_i) X_v+\,\epsilon_i} (MGWR)
#'
#' \eqn{y=\lambda Wy+\beta_c X_c+\,\epsilon_i} (MGWR-SAR(0,k,0))
#'
#' \eqn{y=\lambda Wy+\beta_v(u_i,v_i)X_v+\,\epsilon_i} (MGWR-SAR(0,0,k))
#'
#' \eqn{y=\lambda Wy+\beta_c X_c+\beta_v(u_i,v_i)X_v+\,\epsilon_i} (MGWR-SAR(0,k_c,k_v))
#'
#' \eqn{y=\lambda(u_i,v_i) Wy+\beta_c X_c+\,\epsilon_i} (MGWR-SAR(1,k,0))
#'
#' \eqn{y=\lambda(u_i,v_i)Wy+\beta_v(u_i,v_i)X_v+\,\epsilon_i} (MGWR-SAR(1,0,k))
#'
#' \eqn{y=\lambda(u_i,v_i)Wy+\beta_cX_c+\beta_v(u_i,v_i)X_v+\,\epsilon_i} (MGWR-SAR(1,k_c,k_v))
#'
#' When model imply spatial autocorrelation, a row normalized spatial weight matrix must be provided.
#'2SLS and Best 2SLS method can be used.
#' When model imply local regression, a bandwidth and a kernel type must be provided. Optimal bandwidth can be estimated
#' using bandwidths_mgwrsar function.
#' When model imply mixed local regression, the names of stationary covariates must be provided.
#'
#' #' In addition to the ability of considering spatial autocorrelation in GWR/MGWR like models,
#' MGWRSAR function introduces several useful technics for estimating local regression with space coordinates:
#' \itemize{
#' \item{it uses RCCP and RCCPeigen code that speed up computation and allows parallel computing via doMC package;}
#' \item{it allows to drop out variables with not enough local variance in local regression, which allows to consider dummies in GWR/MGWR framework without trouble.}
#' \item{it allows to drop out local outliers in local regression.}
#' \item{it allows to consider additional variable for kernel, including  time (asymetric kernel) and categorical variables (see Li and Racine 2010). Experimental version.}
#' }
#'
#' @references
#'
#' Geniaux, G. and Martinetti, D. (2017). A new method for dealing simultaneously with spatial autocorrelation and spatial heterogeneity in regression models. Regional Science and Urban Economics. (https://doi.org/10.1016/j.regsciurbeco.2017.04.001)
#'
#' McMillen, D. and Soppelsa, M. E. (2015). A conditionally parametric probit model of
#' microdata land use in chicago. Journal of Regional Science, 55(3):391-415.
#'
#' Loader, C. (1999). Local regression and likelihood, volume 47. springer New York.
#'
#' Franke, R. and Nielson, G. (1980). Smooth interpolation of large sets of scattered data.
#' International journal for numerical methods in engineering, 15(11):1691-1704.
#' @seealso  bandwidths_mgwrsar, summary_mgwrsar, plot_mgwrsar, predict_mgwrsar, kernel_matW
#' @examples
#' \donttest{
#'  library(mgwrsar)
#'  ## loading data example
#'  data(mydata)
#'  coord=as.matrix(mydata[,c("x_lat","y_lon")])
#'  ## Creating a spatial weight matrix (sparce dgCMatrix)
#'  ## of 4 nearest neighbors with 0 in diagonal
#'  W=kernel_matW(H=4,kernels='rectangle',coord_i=coord,NN=4,adaptive=TRUE,
#'  diagnull=TRUE,rowNorm=TRUE)
#'  mgwrsar_0_kc_kv<-MGWRSAR(formula = 'Y_mgwrsar_0_kc_kv~X1+X2+X3', data = mydata,
#'  coord=coord, fixed_vars='X2',kernels=c('gauss'),H=20, Model = 'MGWRSAR_0_kc_kv',
#'  control=list(SE=FALSE,adaptive=TRUE,W=W))
#'  summary_mgwrsar(mgwrsar_0_kc_kv)
#' }
MGWRSAR <- function(formula, data, coord, fixed_vars = NULL, kernels, H,Model = "GWR", control = list()){
  set.seed(123)
 while(sum(duplicated(coord))>0) {
    coord<-jitter(coord,0.001)
    warning('coords have been jittered because there is some duplicated location.')
 }
  colnames(coord)<-c('x_lat','y_lon')
    ptm<-proc.time()
    mycall <- match.call()
    n <-nrow(data)
    assign_control(control,n)
    gwrenv=environment()
    prep_var(gwrenv)
    if (Model == "OLS") {
      model<-list()
      lml <- lm.fit(as.matrix(X), as.matrix(Y))
      model$Betac <- lml$coefficients
      if (SE) {
        rss <- sum(lml$residuals^2)
        rdf <- length(Y) - ncol(X)
        resvar <- rss/rdf
        R <- chol2inv(lml$qr$qr)
        diagR=diag(R)
        model$se =  sqrt(diagR * resvar)
      }
      model$edf = rdf
      model$XC = X
      model$XV = NULL
    }
    else if (Model == "SAR") {
      model<-list()
      if (Method %in% c("B2SLS", "2SLS")) {
        keep = which(!is.na(coefficients(lm.fit(X, Y))))
        mymodel <- mod(as.matrix(Y), as.matrix(X[, keep]),
                       W, as.matrix(X), as.matrix(Y), rep(1, n), "L0",
                       (Method == "B2SLS"), FALSE, SE)
        Betac = mymodel$Betav
        if (SE)
          se = mymodel$se
        else se = NULL
        model$Betac[c(keep, ncol(X) + 1)] <- Betac
        if (SE) {
          model$se <- rep(0, ncol(X) + 1)
          model$se[c(keep, ncol(X) + 1)] <- se
        }
      }
      model$edf = n - ncol(X) - 1
      model$XC = as.matrix(cbind(X, as.matrix(W %*% Y)))
      model$XV = NULL
      model$Betav = NULL
      model$Y = Y
    }
    else if (Model %in%  c("GWR","GWR_glmboost",'GWR_gamboost_linearized',"MGWRSAR_1_0_kv","GWR_multiscale") ){
      if(Model=="MGWRSAR_1_0_kv") Wx=W else Wx=NULL
      model = GWR(Y=Y,XV=XV,ALL_X=X,S=S,H=H,NN=NN,W=Wx, kernels=MykernelS,adaptive=adaptive, Type = Type,SE=SE, isgcv=isgcv,remove_local_outlier=remove_local_outlier,outv=outv,TP=TP,doMC=doMC,ncore=ncore,dists=dists,indexG=indexG,Wd=Wd,Model=Model,S_out=S_out,mstop=mstop,nu=nu)
      model$Betac = NULL
      model$XC = NULL
      if(Model=="MGWRSAR_1_0_kv") {
        if(is.null(new_data)) model$XV = cbind(XV, as.numeric(W %*% Y)) else model$XV = new_XV
        model$edf = model$edf - 1
      } else if(!S_out) model$XV = XV else model$XV = XV[TP,]
    } else {
      if(Model=="MGWR") Wx=NULL else Wx=W
      model<- MGWR(Y=Y,XC=XC,XV=XV,S=S,H=H,NN=NN, kernels=kernels,adaptive=adaptive, Type = Type,SE=SE, isgcv=isgcv,W=W,remove_local_outlier=remove_local_outlier,outv=outv,TP=TP,Model=Model,dists=dists,indexG=indexG,Wd=Wd,S_out=S_out,doMC=doMC,ncore=ncore)
      model$Y = Y
      model$XC=XC
      if(Model =="MGWRSAR_1_kc_kv") {
        if(is.null(new_data)) model$XV = cbind(XV, as.numeric(W %*% Y)) else {
          model$XV = new_XV
          model$XC = new_XC
        }
        model$edf = model$edf - 1
      } else if (Model =="MGWRSAR_1_kc_0") {
        if(is.null(new_data)) model$XV = as.numeric(W %*% Y) else {
          model$XV = NULL
          model$XC = new_XC
        }
      } else model$XV = XV
      if(Model=="MGWRSAR_0_kc_kv") {
        if(is.null(new_data)) model$XC = cbind(XC, as.numeric(W %*% Y)) else {
          model$XV = new_XV
          model$XC = new_XC
        }
        model$edf = model$edf - 1
      } else if(Model=="MGWRSAR_0_0_kv") {
        if(is.null(new_data)) model$XC = as.numeric(W %*% Y) else {
          model$XV = new_XV
          model$XC = NULL
        }
      }
    }
    term1 = 0
    term2 = 0
    if (!is.null(model$Betav) & is.null(new_data)) term1 <- rowSums(model$XV * model$Betav)
    if (!is.null(model$Betac) & is.null(new_data)) term2 <- model$XC %*% as.matrix(model$Betac)
    if(is.null(new_data)) {residuals <- Y - term1 - term2
    fit = as.numeric(term1 + term2)
    }

    if (!is.null(new_data)){
      if(!(Model %in% c('GWR','MGWR'))){
      if(Model %in% c('MGWRSAR_1_0_kv','MGWRSAR_1_kc_kv'))  {
        lambda_pred=model$Betav[,ncol(model$Betav)]
        if(Model =='MGWRSAR_1_0_kv') beta_pred=model$Betav[,-ncol(model$Betav)] else beta_pred=cbind(matrix(model$Betac,nrow=nrow(model$Betav),ncol=length(model$Betac),byrow=TRUE),model$Betav[,-ncol(model$Betav)])
        } else if(Model %in% c('MGWRSAR_0_0_kv','MGWRSAR_0_kc_kv')) {
          lambda_pred=model$Betac[length(model$Betac)]
          if(Model =='MGWRSAR_0_0_kv')  beta_pred=model$Betav else if(Model=='MGWRSAR_0_kc_kv') beta_pred=cbind(matrix(model$Betac[-length(model$Betac)],nrow=nrow(model$Betav),ncol=length(model$Betac)-1,byrow=TRUE),model$Betav)
        } else if(Model =='MGWRSAR_1_kc_0') {
          lambda_pred=model$Betav[,ncol(model$Betav)]
          beta_pred=matrix(model$Betac,nrow=nrow(model$Betav),ncol=length(model$Betac),byrow=TRUE)
        } else if(Model =='SAR') {
          lambda_pred=model$Betac[length(model$Betac)]
          beta_pred=model$Betac[-length(model$Betac)]
        }
        pred=list(lambda_pred=lambda_pred,beta_pred=beta_pred)
        } else  {
        term1 <- rowSums(new_XV * model$Betav)
        if(Model=='MGWR') term2 <- new_XC %*% as.matrix(model$Betac)
        pred = as.numeric(term1 + term2)
      }
    }

    try(colnames(model$Betav) <- names_betav, silent = TRUE)
    try(colnames(model$SEV) <- names_betav, silent = TRUE)
    try(names(model$Betac) <- names_betac, silent = TRUE)
    try(names(model$se) <- names_betac, silent = TRUE)
    if(!is.null(new_data)) {z <- list(Betav = model$Betav, Betac = model$Betac,pred = pred,XC = model$XC, XV = model$XV)} else {
      z <- list(Betav = model$Betav, Betac = model$Betac, sev = model$SEV,se = model$se, Model = Model, Y = Y, XC = model$XC, XV = model$XV,X = X, W = W, isgcv = isgcv, edf = model$edf, tS = model$tS,formula = formula, data = data, Method = Method, coord = coord,H = H, fixed_vars = fixed_vars, kernels = kernels,adaptive=adaptive, fit = fit,pred=pred,residuals = residuals, RMSE=sqrt(mean((residuals[TP])^2)),RMSEn=sqrt(mean((residuals)^2)),SSR = sum(residuals^2), Type = Type,S = S,NN=NN,TP=TP, doMC=doMC,ncore=ncore,mycall = mycall,ctime=(proc.time()-ptm)[3])
    class(z) <- "mgwrsar"}
    invisible(z)
  }

Try the mgwrsar package in your browser

Any scripts or data that you put into this service are public.

mgwrsar documentation built on April 17, 2023, 9:09 a.m.