R/gw.glm.R

In lctools: Local and Geographically Weighted Spatial Statistics Tools

gw.glm<-function(formula, family, dframe, bw, kernel, coords, verbose = TRUE)
  {  
    Obs <- nrow(dframe)

    if (verbose) message("\nNumber of Observations: ", Obs)
    
    if(kernel == 'adaptive')
    { 
      Ne <- bw
      if (verbose) message("\nKernel: Adaptive\nNeighbours: ", Ne)
    } 
    else 
    { 
      if(kernel == 'fixed')
      {
        if (verbose) message("\nKernel: Fixed\nBandwidth: ", bw)
      }
    }
    
    Gl.Model <- eval(substitute(glm(formula, family=family, data = dframe)))
    
    RNames <- names(Gl.Model$coefficients)
    PValueNames <- paste("P",RNames, sep="_")

    ModelVarNo <- length(Gl.Model$coefficients)
    
    if (verbose) {
      message("\nNumber of Variables: ", ModelVarNo-1)
      message("\n--------------- Global Model Summary ---------------\n")
      print(summary(Gl.Model))
      }
    
    Gl.Model.res<-Gl.Model$y-Gl.Model$fitted.values
    Gl.RSS<-sum(Gl.Model.res^2)
    mean.y<-mean(Gl.Model$y)
    TSS<-sum((Gl.Model$y-mean.y)^2)
    
    Gl.r2<-1-(Gl.RSS/TSS)
    Gl.r2.adj<-1-((Gl.RSS/(Obs-ModelVarNo-1))/(TSS/(Obs-1)))
    
    if (verbose) {
      message("\nResidual Sum of Squares: ", Gl.RSS)
      message("\nR-squared: ", Gl.r2)
      message("\nAdjusted R-squared: ", Gl.r2.adj)
    }
    
    DistanceT <- dist(coords)
    Dij <- as.matrix(DistanceT)

    GGLM_LEst <- as.data.frame(setNames(replicate(ModelVarNo,numeric(0), simplify = F), RNames[1:ModelVarNo]))
    GGLM_LPvalues <- as.data.frame(setNames(replicate(ModelVarNo,numeric(0), simplify = F), PValueNames[1:ModelVarNo]))
    GGLM_GofFit <- data.frame(y=numeric(0), GLM_yfit=numeric(0), GLM_Res=numeric(0), 
                            GLM_AIC=numeric(0), GLM_Deviance=numeric(0),GLM_DevExpl=numeric(0),
                            GLM_Overdispersion=numeric(0),GLM_OverDis1=numeric(0))
  
  for(m in 1:Obs){
    #Get the data
    DNeighbour <- Dij[,m]
    DataSet <- data.frame(dframe, DNeighbour=DNeighbour)

    #Sort by distance
    DataSetSorted <- DataSet[order(DataSet$DNeighbour),]

    if(kernel == 'adaptive')
      { 
        #Keep Nearest Neighbours
        SubSet <- DataSetSorted[1:Ne,]
        Kernel_H <- max(SubSet$DNeighbour)
      } 
    else 
      { 
        if(kernel == 'fixed')
          {
             SubSet <- subset(DataSetSorted, DNeighbour <= bw)
             Kernel_H <- bw
          }
      }
    
    #Bi-square weights
    Wts<-(1-(SubSet$DNeighbour/Kernel_H)^2)^2
    
    #Calculate WGLM
    Lcl.Model<-eval(substitute(glm(formula, family=family, data = SubSet, weights=Wts)))
    
    #Store in table
    GGLM_LEst[m,]<-Lcl.Model$coefficients
    GGLM_LPvalues[m,]<-summary(Lcl.Model)$coefficients[,4]
    
    GGLM_GofFit[m,1]<-Gl.Model$y[m]
    GGLM_GofFit[m,2]<-exp(sum(Gl.Model$model[m,2:ModelVarNo] * Lcl.Model$coefficients[2:ModelVarNo]) + sum(Lcl.Model$coefficients[1]))
   
    GGLM_GofFit[m,4]<-AIC(Lcl.Model)
    GGLM_GofFit[m,5]<-deviance(Lcl.Model)
    GGLM_GofFit[m,6]<-100*((Lcl.Model$null.deviance-Lcl.Model$deviance)/Lcl.Model$null.deviance)
    GGLM_GofFit[m,7]<-Lcl.Model$deviance/Lcl.Model$df.residual
    GGLM_GofFit[m,8]<-var(Lcl.Model$model[,1])/mean(Lcl.Model$model[,1])
  
  }
   
    Adj.f.l <- sum(GGLM_GofFit[,1])/sum(GGLM_GofFit[,2])
    GGLM_GofFit[,3] <- GGLM_GofFit[,1] - Adj.f.l*GGLM_GofFit[,2]
    
  gw.glm.out<-list(GGLM_LEst=GGLM_LEst, GGLM_LPvalues=GGLM_LPvalues, GGLM_GofFit=GGLM_GofFit)
  if (verbose) {
    message("\n--------------- Local Model Summary ---------------\n")
    message("\nResiduals:\n")
    print(summary(gw.glm.out$GGLM_GofFit$GLM_Res))
  }
  
  t1 <- data.frame(Min = apply(gw.glm.out$GGLM_LEst, 2, min), Max = apply(gw.glm.out$GGLM_LEst, 2, max), 
            Mean = apply(gw.glm.out$GGLM_LEst, 2, mean), StD = apply(gw.glm.out$GGLM_LEst, 2, sd))
  
  if (verbose) {
    message("\nCoefficients:\n")
    print(t1)
  }

  l.RSS<-sum(gw.glm.out$GGLM_GofFit$GLM_Res^2)
  
  l.r2<-1-(l.RSS/TSS)
  l.r2.adj<-1-((l.RSS/(Obs-ModelVarNo-1))/(TSS/(Obs-1)))
  
  if (verbose) {
    message("\nResidual Sum of Squares: ", l.RSS)
    message("\nR-squared: ", l.r2)
    message("\nAdjusted R-squared: ", l.r2.adj)
  }
  
  return(gw.glm.out)
}