Final Pull Request for adding Poisson MGWR

mgwr/gwr.py

@@ -134,7 +134,7 @@ class GWR(GLM):
     spherical     : boolean
                     True for shperical coordinates (long-lat),
                     False for projected coordinates (defalut).
-                    
+
     hat_matrix    : boolean
                     True to store full n by n hat matrix,
                     False to not store full hat matrix to minimize memory footprint (defalut).
@@ -516,10 +516,10 @@ class GWRResults(GLMResults):
 
     R2                  : float
                           R-squared for the entire model (1- RSS/TSS)
-                          
+
     adj_R2              : float
                           adjusted R-squared for the entire model
-                          
+
     aic                 : float
                           Akaike information criterion
 
@@ -575,11 +575,11 @@ class GWRResults(GLMResults):
     pDev                : float
                           local percent of deviation accounted for; analogous to
                           r-squared for GLM's
-                          
+
     D2                  : float
                           percent deviance explained for GLM, equivaleng to R2 for
                           Gaussian.
-                          
+
     adj_D2              : float
                           adjusted percent deviance explained, equivaleng to adjusted
                           R2 for Gaussian.
@@ -1438,8 +1438,8 @@ class MGWR(GWR):
 
     """
 
-    def __init__(self, coords, y, X, selector, sigma2_v1=True,
-                 kernel='bisquare', fixed=False, constant=True,
+    def __init__(self, coords, y, X, selector,family=Gaussian(), offset=None,
+                 sigma2_v1=True, kernel='bisquare', fixed=False, constant=True,
                  spherical=False, hat_matrix=False):
         """
         Initialize class
@@ -1448,17 +1448,20 @@ def __init__(self, coords, y, X, selector, sigma2_v1=True,
         self.bws = self.selector.bw[0]  #final set of bandwidth
         self.bws_history = selector.bw[1]  #bws history in backfitting
         self.bw_init = self.selector.bw_init  #initialization bandiwdth
-        self.family = Gaussian(
-        )  # manually set since we only support Gassian MGWR for now
+        self.family = family
+        if offset is None:
+            self.offset = np.ones((len(y), 1))
+        else:
+            self.offset = offset * 1.0
         GWR.__init__(self, coords, y, X, self.bw_init, family=self.family,
-                     sigma2_v1=sigma2_v1, kernel=kernel, fixed=fixed,
-                     constant=constant, spherical=spherical,
+                     offset=self.offset, sigma2_v1=sigma2_v1, kernel=kernel,
+                     fixed=fixed, constant=constant, spherical=spherical,
                      hat_matrix=hat_matrix)
         self.selector = selector
         self.sigma2_v1 = sigma2_v1
         self.points = None
         self.P = None
-        self.offset = None
+        self.family = family
         self.exog_resid = None
         self.exog_scale = None
         self_fit_params = None
@@ -1483,7 +1486,14 @@ def _chunk_compute_R(self, chunk_id=0):
 
         for i in range(n):
             wi = self._build_wi(i, self.bw_init).reshape(-1, 1)
-            xT = (self.X * wi).T
+            if isinstance(self.family, Poisson):
+                wi = wi.reshape(-1,1)
+                rslt = iwls(self.y, self.X, self.family, self.offset, None,
+                            wi=wi)
+                w = rslt[3]
+                xT = (self.X * w).T
+            else:
+                xT = (self.X * wi).T
             P = np.linalg.solve(xT.dot(self.X), xT).dot(init_pR).T
             pR[i, :, :] = P * self.X[i]
 
@@ -1502,8 +1512,16 @@ def _chunk_compute_R(self, chunk_id=0):
                     for i in range(len(chunk_index_Aj)):
                         index = chunk_index_Aj[i]
                         wi = self._build_wi(index, self.bws_history[iter_i, j])
-                        xw = Xj * wi
-                        pAj[i, :] = Xj[index] / np.sum(xw * Xj) * xw
+                        if isinstance(self.family, Poisson):
+                            Xj = Xj.reshape(-1,1)
+                            wi = wi.reshape(-1,1)
+                            rslt = iwls(self.y, Xj, self.family, self.offset, None,
+                                        wi=wi)
+                            w = rslt[3]
+                            xw = Xj * w
+                        else:
+                            xw = Xj * wi
+                        pAj[i, :] = (Xj[index] / np.sum(xw * Xj) * xw).reshape(-1)
                     pR[chunk_index_Aj, :, j] = pAj.dot(pRj_old)
                 err = pRj_old - pR[:, :, j]
 
@@ -1520,21 +1538,25 @@ def _chunk_compute_R(self, chunk_id=0):
     def fit(self, n_chunks=1, pool=None):
         """
         Compute MGWR inference by chunk to reduce memory footprint.
-        
+
         Parameters
         ----------
 
         n_chunks      : integer, optional
-                        A number of chunks parameter to reduce memory usage. 
+                        A number of chunks parameter to reduce memory usage.
                         e.g. n_chunks=2 should reduce overall memory usage by 2.
         pool          : A multiprocessing Pool object to enable parallel fitting; default is None.
-                        
+
         Returns
         -------
                       : MGWRResults
         """
         params = self.selector.params
-        predy = np.sum(self.X * params, axis=1).reshape(-1, 1)
+
+        if isinstance(self.family, Poisson):
+            predy = self.offset*(np.exp(np.sum(self.X * params, axis=1).reshape(-1, 1)))
+        else:
+            predy = np.sum(self.X * params, axis=1).reshape(-1, 1)
 
         try:
             from tqdm.autonotebook import tqdm  #progress bar
@@ -1692,7 +1714,7 @@ class MGWRResults(GWRResults):
 
     R2                  : float
                           R-squared for the entire model (1- RSS/TSS)
-                          
+
     adj_R2              : float
                           adjusted R-squared for the entire model
 

mgwr/search.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 from copy import deepcopy
+from spglm.family import Gaussian, Binomial, Poisson
 
 
 def golden_section(a, c, delta, function, tol, max_iter, int_score=False,
@@ -164,7 +165,7 @@ def equal_interval(l_bound, u_bound, interval, function, int_score=False,
     return opt_val, opt_score, output
 
 
-def multi_bw(init, y, X, n, k, family, tol, max_iter, rss_score, gwr_func,
+def multi_bw(init, y, X, n, k, family, offset, tol, max_iter, rss_score, gwr_func,
              bw_func, sel_func, multi_bw_min, multi_bw_max, bws_same_times,
              verbose=False):
     """
@@ -180,7 +181,10 @@ def multi_bw(init, y, X, n, k, family, tol, max_iter, rss_score, gwr_func,
     err = optim_model.resid_response.reshape((-1, 1))
     param = optim_model.params
 
-    XB = np.multiply(param, X)
+    if isinstance(family, Poisson):
+        XB = offset*np.exp(np.multiply(param, X))
+    else:
+        XB = np.multiply(param, X)
     if rss_score:
         rss = np.sum((err)**2)
     iters = 0
@@ -230,7 +234,10 @@ def tqdm(x, desc=''):  #otherwise, just passthrough the range
         XB = new_XB
 
         if rss_score:
-            predy = np.sum(np.multiply(params, X), axis=1).reshape((-1, 1))
+            if isinstance(family, Poisson):
+                predy = offset*(np.exp(np.sum(X * params, axis=1).reshape(-1, 1)))
+            else:
+                predy = np.sum(np.multiply(params, X), axis=1).reshape((-1, 1))
             new_rss = np.sum((y - predy)**2)
             score = np.abs((new_rss - rss) / new_rss)
             rss = new_rss

mgwr/sel_bw.py

@@ -136,7 +136,7 @@ class Sel_BW(object):
     >>> pov = np.array(data.by_col('PctPov')).reshape((-1,1))
     >>> african_amer = np.array(data.by_col('PctBlack')).reshape((-1,1))
     >>> X = np.hstack([rural, pov, african_amer])
-    
+
     Golden section search AICc - adaptive bisquare
 
     >>> bw = Sel_BW(coords, y, X).search(criterion='AICc')
@@ -213,7 +213,7 @@ def search(self, search_method='golden_section', criterion='AICc',
                          min value used in bandwidth search
         bw_max         : float
                          max value used in bandwidth search
-        multi_bw_min   : list 
+        multi_bw_min   : list
                          min values used for each covariate in mgwr bandwidth search.
                          Must be either a single value or have one value for
                          each covariate including the intercept
@@ -391,10 +391,10 @@ def sel_func(bw_func, bw_min=None, bw_max=None):
                 bw_min=bw_min, bw_max=bw_max, interval=interval, tol=tol,
                 max_iter=max_iter, pool=self.pool, verbose=False)
 
-        self.bw = multi_bw(self.init_multi, y, X, n, k, family, self.tol_multi,
-                           self.max_iter_multi, self.rss_score, gwr_func,
-                           bw_func, sel_func, multi_bw_min, multi_bw_max,
-                           bws_same_times, verbose=self.verbose)
+        self.bw = multi_bw(self.init_multi, y, X, n, k, family, offset,
+                           self.tol_multi, self.max_iter_multi, self.rss_score,
+                           gwr_func, bw_func, sel_func, multi_bw_min,
+                           multi_bw_max, bws_same_times, verbose=self.verbose)
 
     def _init_section(self, X_glob, X_loc, coords, constant):
         if len(X_glob) > 0: