-
Notifications
You must be signed in to change notification settings - Fork 121
/
gwr.py
executable file
·2195 lines (1765 loc) · 74.6 KB
/
gwr.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Main GWR classes
__author__ = "Taylor Oshan Tayoshan@gmail.com"
import copy
import os
from typing import Optional
import warnings
import numpy as np
import numpy.linalg as la
from scipy.stats import t
from scipy.special import factorial
from itertools import combinations as combo
from spglm.family import Gaussian, Binomial, Poisson
from spglm.glm import GLM, GLMResults
from spglm.iwls import iwls, _compute_betas_gwr
from spglm.utils import cache_readonly
from joblib import Parallel, delayed
from .diagnostics import get_AIC, get_AICc, get_BIC, corr
from .kernels import *
from .summary import *
class GWR(GLM):
"""
Geographically weighted regression. Can currently estimate Gaussian,
Poisson, and logistic models(built on a GLM framework). GWR object prepares
model input. Fit method performs estimation and returns a GWRResults object.
Parameters
----------
coords : array-like
n*2, collection of n sets of (x,y) coordinates of
observatons; also used as calibration locations is
'points' is set to None
y : array
n*1, dependent variable
X : array
n*k, independent variable, exlcuding the constant
bw : scalar
bandwidth value consisting of either a distance or N
nearest neighbors; user specified or obtained using
Sel_BW
family : family object
underlying probability model; provides
distribution-specific calculations
offset : array
n*1, the offset variable at the ith location. For Poisson model
this term is often the size of the population at risk or
the expected size of the outcome in spatial epidemiology
Default is None where Ni becomes 1.0 for all locations;
only for Poisson models
sigma2_v1 : boolean
specify form of corrected denominator of sigma squared to use for
model diagnostics; Acceptable options are:
'True': n-tr(S) (defualt)
'False': n-2(tr(S)+tr(S'S))
kernel : string
type of kernel function used to weight observations;
available options:
'gaussian'
'bisquare'
'exponential'
fixed : boolean
True for distance based kernel function and False for
adaptive (nearest neighbor) kernel function (default)
constant : boolean
True to include intercept (default) in model and False to exclude
intercept.
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).
name_x : list of strings
Names of independent variables for use in output
n_jobs : integer
The number of jobs (default -1) to run in parallel. -1 means using all processors.
Attributes
----------
coords : array-like
n*2, collection of n sets of (x,y) coordinates used for
calibration locations
y : array
n*1, dependent variable
X : array
n*k, independent variable, exlcuding the constant
bw : scalar
bandwidth value consisting of either a distance or N
nearest neighbors; user specified or obtained using
Sel_BW
family : family object
underlying probability model; provides
distribution-specific calculations
offset : array
n*1, the offset variable at the ith location. For Poisson model
this term is often the size of the population at risk or
the expected size of the outcome in spatial epidemiology
Default is None where Ni becomes 1.0 for all locations
sigma2_v1 : boolean
specify form of corrected denominator of sigma squared to use for
model diagnostics; Acceptable options are:
'True': n-tr(S) (defualt)
'False': n-2(tr(S)+tr(S'S))
kernel : string
type of kernel function used to weight observations;
available options:
'gaussian'
'bisquare'
'exponential'
fixed : boolean
True for distance based kernel function and False for
adaptive (nearest neighbor) kernel function (default)
constant : boolean
True to include intercept (default) in model and False to exclude
intercept
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).
n : integer
number of observations
k : integer
number of independent variables
mean_y : float
mean of y
std_y : float
standard deviation of y
fit_params : dict
parameters passed into fit method to define estimation
routine
points : array-like
n*2, collection of n sets of (x,y) coordinates used for
calibration locations instead of all observations;
defaults to None unles specified in predict method
P : array
n*k, independent variables used to make prediction;
exlcuding the constant; default to None unless specified
in predict method
exog_scale : scalar
estimated scale using sampled locations; defualt is None
unless specified in predict method
exog_resid : array-like
estimated residuals using sampled locations; defualt is None
unless specified in predict method
Examples
--------
#basic model calibration
>>> import libpysal as ps
>>> from mgwr.gwr import GWR
>>> data = ps.io.open(ps.examples.get_path('GData_utm.csv'))
>>> coords = list(zip(data.by_col('X'), data.by_col('Y')))
>>> y = np.array(data.by_col('PctBach')).reshape((-1,1))
>>> rural = np.array(data.by_col('PctRural')).reshape((-1,1))
>>> 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])
>>> model = GWR(coords, y, X, bw=90.000, fixed=False, kernel='bisquare')
>>> results = model.fit()
>>> print(results.params.shape)
(159, 4)
#predict at unsampled locations
>>> index = np.arange(len(y))
>>> test = index[-10:]
>>> X_test = X[test]
>>> coords_test = np.array(coords)[test]
>>> model = GWR(coords, y, X, bw=94, fixed=False, kernel='bisquare')
>>> results = model.predict(coords_test, X_test)
>>> print(results.params.shape)
(10, 4)
"""
def __init__(self, coords, y, X, bw, family=Gaussian(), offset=None,
sigma2_v1=True, kernel='bisquare', fixed=False, constant=True,
spherical=False, hat_matrix=False, name_x=None,n_jobs=-1):
"""
Initialize class
"""
GLM.__init__(self, y, X, family, constant=constant)
self.constant = constant
self.sigma2_v1 = sigma2_v1
self.coords = np.array(coords)
self.bw = bw
self.kernel = kernel
self.fixed = fixed
if offset is None:
self.offset = np.ones((self.n, 1))
else:
self.offset = offset * 1.0
self.fit_params = {}
self.points = None
self.exog_scale = None
self.exog_resid = None
self.P = None
self.spherical = spherical
self.hat_matrix = hat_matrix
self.name_x = name_x
self.n_jobs = n_jobs
def _build_wi(self, i, bw):
if bw == np.inf:
wi = np.ones((self.n))
return wi
try:
wi = Kernel(i, self.coords, bw, fixed=self.fixed,
function=self.kernel, points=self.points,
spherical=self.spherical).kernel
except BaseException:
raise # TypeError('Unsupported kernel function ', kernel)
return wi
def _local_fit(self, i):
"""
Local fitting at location i.
"""
wi = self._build_wi(i, self.bw).reshape(-1, 1) #local spatial weights
if isinstance(self.family, Gaussian):
betas, inv_xtx_xt = _compute_betas_gwr(self.y, self.X, wi)
predy = np.dot(self.X[i], betas)[0]
resid = self.y[i] - predy
influ = np.dot(self.X[i], inv_xtx_xt[:, i])
w = 1
elif isinstance(self.family, (Poisson, Binomial)):
rslt = iwls(self.y, self.X, self.family, self.offset, None,
self.fit_params['ini_params'], self.fit_params['tol'],
self.fit_params['max_iter'], wi=wi)
inv_xtx_xt = rslt[5]
w = rslt[3][i][0]
influ = np.dot(self.X[i], inv_xtx_xt[:, i]) * w
predy = rslt[1][i]
resid = self.y[i] - predy
betas = rslt[0]
if self.fit_params['lite']:
return influ, resid, predy, betas.reshape(-1)
else:
Si = np.dot(self.X[i], inv_xtx_xt).reshape(-1)
tr_STS_i = np.sum(Si * Si * w * w)
CCT = np.diag(np.dot(inv_xtx_xt, inv_xtx_xt.T)).reshape(-1)
if not self.hat_matrix:
Si = None
return influ, resid, predy, betas.reshape(-1), w, Si, tr_STS_i, CCT
def fit(self, ini_params=None, tol=1.0e-5, max_iter=20, solve='iwls',
lite=False,pool=None):
"""
Method that fits a model with a particular estimation routine.
Parameters
----------
ini_betas : array, optional
k*1, initial coefficient values, including constant.
Default is None, which calculates initial values during
estimation.
tol: float, optional
Tolerence for estimation convergence.
Default is 1.0e-5.
max_iter : integer, optional
Maximum number of iterations if convergence not
achieved. Default is 20.
solve : string, optional
Technique to solve MLE equations.
Default is 'iwls', meaning iteratively (
re)weighted least squares.
lite : bool, optional
Whether to estimate a lightweight GWR that
computes the minimum diagnostics needed for
bandwidth selection (could speed up
bandwidth selection for GWR) or to estimate
a full GWR. Default is False.
pool : None, deprecated and not used.
Returns
-------
:
If lite=False, return a GWRResult
instance; otherwise, return a GWRResultLite
instance.
"""
self.fit_params['ini_params'] = ini_params
self.fit_params['tol'] = tol
self.fit_params['max_iter'] = max_iter
self.fit_params['solve'] = solve
self.fit_params['lite'] = lite
if pool:
warnings.warn("The pool parameter is no longer used and will have no effect; parallelization is default and implemented using joblib instead.", RuntimeWarning, stacklevel=2)
if solve.lower() == 'iwls':
if self.points is None:
m = self.y.shape[0]
else:
m = self.points.shape[0]
rslt = Parallel(n_jobs=self.n_jobs)(delayed(self._local_fit)(i) for i in range(m))
rslt_list = list(zip(*rslt))
influ = np.array(rslt_list[0]).reshape(-1, 1)
resid = np.array(rslt_list[1]).reshape(-1, 1)
params = np.array(rslt_list[3])
if lite:
return GWRResultsLite(self, resid, influ, params)
else:
predy = np.array(rslt_list[2]).reshape(-1, 1)
w = np.array(rslt_list[-4]).reshape(-1, 1)
if self.hat_matrix:
S = np.array(rslt_list[-3])
else:
S = None
tr_STS = np.sum(np.array(rslt_list[-2]))
CCT = np.array(rslt_list[-1])
return GWRResults(self, params, predy, S, CCT, influ, tr_STS,
w, self.name_x)
def predict(self, points, P, exog_scale=None, exog_resid=None,
fit_params={}):
"""
Method that predicts values of the dependent variable at un-sampled
locations
Parameters
----------
points : array-like
n*2, collection of n sets of (x,y) coordinates used for
calibration prediction locations
P : array
n*k, independent variables used to make prediction;
exlcuding the constant
exog_scale : scalar
estimated scale using sampled locations; defualt is None
which estimates a model using points from "coords"
exog_resid : array-like
estimated residuals using sampled locations; defualt is None
which estimates a model using points from "coords"; if
given it must be n*1 where n is the length of coords
fit_params : dict
key-value pairs of parameters that will be passed into fit
method to define estimation routine; see fit method for more details
"""
if (exog_scale is None) & (exog_resid is None):
train_gwr = self.fit(**fit_params)
self.exog_scale = train_gwr.scale
self.exog_resid = train_gwr.resid_response
elif (exog_scale is not None) & (exog_resid is not None):
self.exog_scale = exog_scale
self.exog_resid = exog_resid
else:
raise InputError('exog_scale and exog_resid must both either be'
'None or specified')
self.points = points
if self.constant:
P = np.hstack([np.ones((len(P), 1)), P])
self.P = P
else:
self.P = P
gwr = self.fit(**fit_params)
return gwr
@cache_readonly
def df_model(self):
return None
@cache_readonly
def df_resid(self):
return None
class GWRResults(GLMResults):
"""
Basic class including common properties for all GWR regression models
Parameters
----------
model : GWR object
pointer to GWR object with estimation parameters
params : array
n*k, estimated coefficients
predy : array
n*1, predicted y values
S : array
n*n, hat matrix
CCT : array
n*k, scaled variance-covariance matrix
w : array
n*1, final weight used for iteratively re-weighted least
sqaures; default is None
name_x : list of strings
Names of independent variables for use in output
Attributes
----------
model : GWR Object
points to GWR object for which parameters have been
estimated
params : array
n*k, parameter estimates
predy : array
n*1, predicted value of y
y : array
n*1, dependent variable
X : array
n*k, independent variable, including constant
family : family object
underlying probability model; provides
distribution-specific calculations
n : integer
number of observations
k : integer
number of independent variables
df_model : integer
model degrees of freedom
df_resid : integer
residual degrees of freedom
offset : array
n*1, the offset variable at the ith location.
For Poisson model this term is often the size of
the population at risk or the expected size of
the outcome in spatial epidemiology; Default is
None where Ni becomes 1.0 for all locations
scale : float
sigma squared used for subsequent computations
w : array
n*1, final weights from iteratively re-weighted least
sqaures routine
resid_response : array
n*1, residuals of the repsonse
resid_ss : scalar
residual sum of sqaures
W : array
n*n; spatial weights for each observation from each
calibration point
S : array
n*n, hat matrix
CCT : array
n*k, scaled variance-covariance matrix
ENP : scalar
effective number of paramters, which depends on
sigma2
tr_S : float
trace of S (hat) matrix
tr_STS : float
trace of STS matrix
y_bar : array
n*1, weighted mean value of y
TSS : array
n*1, geographically weighted total sum of squares
RSS : array
n*1, geographically weighted residual sum of squares
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
aicc : float
corrected Akaike information criterion to account
to account for model complexity (smaller
bandwidths)
bic : float
Bayesian information criterio
localR2 : array
n*1, local R square
sigma2 : float
sigma squared (residual variance) that has been
corrected to account for the ENP
std_res : array
n*1, standardised residuals
bse : array
n*k, standard errors of parameters (betas)
influ : array
n*1, leading diagonal of S matrix
CooksD : array
n*1, Cook's D
tvalues : array
n*k, local t-statistics
adj_alpha : array
3*1, corrected alpha values to account for multiple
hypothesis testing for the 90%, 95%, and 99% confidence
levels; tvalues with an absolute value larger than the
corrected alpha are considered statistically
significant.
deviance : array
n*1, local model deviance for each calibration point
resid_deviance : array
n*1, local sum of residual deviance for each
calibration point
llf : scalar
log-likelihood of the full model; see
pysal.contrib.glm.family for damily-sepcific
log-likelihoods
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.
mu : array
n*, flat one dimensional array of predicted mean
response value from estimator
fit_params : dict
parameters passed into fit method to define estimation
routine
predictions : array
p*1, predicted values generated by calling the GWR
predict method to predict dependent variable at
unsampled points ()
name_x : list of strings
Names of independent variables for use in output
"""
def __init__(self, model, params, predy, S, CCT, influ, tr_STS=None,
w=None, name_x=None):
GLMResults.__init__(self, model, params, predy, w)
self.offset = model.offset
if w is not None:
self.w = w
self.predy = predy
self.S = S
self.tr_STS = tr_STS
self.influ = influ
self.CCT = self.cov_params(CCT, model.exog_scale)
self._cache = {}
self.name_x=name_x
@cache_readonly
def W(self):
W = np.array(
[self.model._build_wi(i, self.model.bw) for i in range(self.n)])
return W
@cache_readonly
def resid_ss(self):
if self.model.points is not None:
raise NotImplementedError('Not available for GWR prediction')
else:
u = self.resid_response.flatten()
return np.dot(u, u.T)
@cache_readonly
def scale(self, scale=None):
if isinstance(self.family, Gaussian):
scale = self.sigma2
else:
scale = 1.0
return scale
def cov_params(self, cov, exog_scale=None):
"""
Returns scaled covariance parameters
Parameters
----------
cov : array
estimated covariance parameters
Returns
-------
Scaled covariance parameters
"""
if exog_scale is not None:
return cov * exog_scale
else:
return cov * self.scale
@cache_readonly
def tr_S(self):
"""
trace of S (hat) matrix
"""
return np.sum(self.influ)
@cache_readonly
def ENP(self):
"""
effective number of parameters
Defaults to tr(s) as defined in :cite:`yu:2019`
but can alternatively be based on 2tr(s) - tr(STS)
and the form depends on the specification of sigma2
"""
if self.model.sigma2_v1:
return self.tr_S
else:
return 2 * self.tr_S - self.tr_STS
@cache_readonly
def y_bar(self):
"""
weighted mean of y
"""
if self.model.points is not None:
n = len(self.model.points)
else:
n = self.n
off = self.offset.reshape((-1, 1))
arr_ybar = np.zeros(shape=(self.n, 1))
for i in range(n):
w_i = np.reshape(self.model._build_wi(i, self.model.bw), (-1, 1))
sum_yw = np.sum(self.y.reshape((-1, 1)) * w_i)
arr_ybar[i] = 1.0 * sum_yw / np.sum(w_i * off)
return arr_ybar
@cache_readonly
def TSS(self):
"""
geographically weighted total sum of squares
Methods: p215, (9.9)
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
"""
if self.model.points is not None:
n = len(self.model.points)
else:
n = self.n
TSS = np.zeros(shape=(n, 1))
for i in range(n):
TSS[i] = np.sum(
np.reshape(self.model._build_wi(i, self.model.bw),
(-1, 1)) * (self.y.reshape(
(-1, 1)) - self.y_bar[i])**2)
return TSS
@cache_readonly
def RSS(self):
"""
geographically weighted residual sum of squares
Methods: p215, (9.10)
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
"""
if self.model.points is not None:
n = len(self.model.points)
resid = self.model.exog_resid.reshape((-1, 1))
else:
n = self.n
resid = self.resid_response.reshape((-1, 1))
RSS = np.zeros(shape=(n, 1))
for i in range(n):
RSS[i] = np.sum(
np.reshape(self.model._build_wi(i, self.model.bw),
(-1, 1)) * resid**2)
return RSS
@cache_readonly
def localR2(self):
"""
local R square
Methods: p215, (9.8)
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
"""
if isinstance(self.family, Gaussian):
return (self.TSS - self.RSS) / self.TSS
else:
raise NotImplementedError('Only applicable to Gaussian')
@cache_readonly
def sigma2(self):
"""
residual variance
if sigma2_v1 is True: only use n-tr(S) in denominator
Methods: p214, (9.6) :cite:`fotheringham_geographically_2002`
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
and as defined in :cite:`yu:2019`
if sigma2_v1 is False (v1v2): use n-2(tr(S)+tr(S'S)) in denominator
Methods: p55 (2.16)-(2.18) :cite:`fotheringham_geographically_2002`
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
"""
if self.model.sigma2_v1:
return (self.resid_ss / (self.n - self.tr_S))
else:
# could be changed to SWSTW - nothing to test against
return self.resid_ss / (self.n - 2.0 * self.tr_S + self.tr_STS)
@cache_readonly
def std_res(self):
"""
standardized residuals
Methods: p215, (9.7) :cite:`fotheringham_geographically_2002`
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
"""
return self.resid_response.reshape(
(-1, 1)) / (np.sqrt(self.scale * (1.0 - self.influ)))
@cache_readonly
def bse(self):
"""
standard errors of Betas
Methods: p215, (2.15) and (2.21) :cite:`fotheringham_geographically_2002`
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
"""
return np.sqrt(self.CCT)
@cache_readonly
def cooksD(self):
"""
Influence: leading diagonal of S Matrix
Methods: p216, (9.11) :cite:`fotheringham_geographically_2002`
Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002).
Geographically weighted regression: the analysis of spatially varying
relationships.
Note: in (9.11), p should be tr(S), that is, the effective number of parameters
"""
return self.std_res**2 * self.influ / (self.tr_S * (1.0 - self.influ))
@cache_readonly
def deviance(self):
off = self.offset.reshape((-1, 1)).T
y = self.y
ybar = self.y_bar
if isinstance(self.family, Gaussian):
raise NotImplementedError(
'deviance not currently used for Gaussian')
elif isinstance(self.family, Poisson):
dev = np.sum(
2.0 * self.W * (y * np.log(y / (ybar * off)) -
(y - ybar * off)), axis=1)
elif isinstance(self.family, Binomial):
dev = self.family.deviance(self.y, self.y_bar, self.W, axis=1)
return dev.reshape((-1, 1))
@cache_readonly
def resid_deviance(self):
if isinstance(self.family, Gaussian):
raise NotImplementedError(
'deviance not currently used for Gaussian')
else:
off = self.offset.reshape((-1, 1)).T
y = self.y
ybar = self.y_bar
global_dev_res = ((self.family.resid_dev(self.y, self.mu))**2)
dev_res = np.repeat(global_dev_res.flatten(), self.n)
dev_res = dev_res.reshape((self.n, self.n))
dev_res = np.sum(dev_res * self.W.T, axis=0)
return dev_res.reshape((-1, 1))
@cache_readonly
def pDev(self):
"""
Local percentage of deviance accounted for. Described in the GWR4
manual. Equivalent to 1 - (deviance/null deviance)
"""
if isinstance(self.family, Gaussian):
raise NotImplementedError('Not implemented for Gaussian')
else:
return 1.0 - (self.resid_deviance / self.deviance)
@cache_readonly
def adj_alpha(self):
"""
Corrected alpha (critical) values to account for multiple testing during hypothesis
testing. Includes corrected value for 90% (.1), 95% (.05), and 99%
(.01) confidence levels. Correction comes from:
:cite:`Silva:2016` : da Silva, A. R., & Fotheringham, A. S. (2015). The Multiple Testing Issue in
Geographically Weighted Regression. Geographical Analysis.
"""
alpha = np.array([.1, .05, .001])
pe = self.ENP
p = self.k
return (alpha * p) / pe
def critical_tval(self, alpha=None):
"""
Utility function to derive the critical t-value based on given alpha
that are needed for hypothesis testing
Parameters
----------
alpha : scalar
critical value to determine which tvalues are
associated with statistically significant parameter
estimates. Default to None in which case the adjusted
alpha value at the 95 percent CI is automatically
used.
Returns
-------
critical : scalar
critical t-val based on alpha
"""
n = self.n
if alpha is not None:
alpha = np.abs(alpha) / 2.0
critical = t.ppf(1 - alpha, n - 1)
else:
alpha = np.abs(self.adj_alpha[1]) / 2.0
critical = t.ppf(1 - alpha, n - 1)
return critical
def filter_tvals(self, critical_t=None, alpha=None):
"""
Utility function to set tvalues with an absolute value smaller than the
absolute value of the alpha (critical) value to 0. If critical_t
is supplied than it is used directly to filter. If alpha is provided
than the critical t value will be derived and used to filter. If neither
are critical_t nor alpha are provided, an adjusted alpha at the 95
percent CI will automatically be used to define the critical t-value and
used to filter. If both critical_t and alpha are supplied then the alpha
value will be ignored.
Parameters
----------
critical_t : scalar
critical t-value to determine whether parameters are
statistically significant
alpha : scalar
alpha value to determine which tvalues are
associated with statistically significant parameter
estimates
Returns
-------
filtered : array
n*k; new set of n tvalues for each of k variables
where absolute tvalues less than the absolute value of
alpha have been set to 0.
"""
n = self.n
if critical_t is not None:
critical = critical_t
else:
critical = self.critical_tval(alpha=alpha)
subset = (self.tvalues < critical) & (self.tvalues > -1.0 * critical)
tvalues = self.tvalues.copy()
tvalues[subset] = 0
return tvalues
@cache_readonly
def df_model(self):
return self.n - self.tr_S
@cache_readonly
def df_resid(self):
return self.n - 2.0 * self.tr_S + self.tr_STS
@cache_readonly
def normalized_cov_params(self):
return None
@cache_readonly
def resid_pearson(self):
return None
@cache_readonly
def resid_working(self):
return None
@cache_readonly
def resid_anscombe(self):
return None
@cache_readonly
def pearson_chi2(self):
return None
@cache_readonly
def llnull(self):