1""" 2Covariance models and estimators for GEE. 3 4Some details for the covariance calculations can be found in the Stata 5docs: 6 7http://www.stata.com/manuals13/xtxtgee.pdf 8""" 9from statsmodels.compat.pandas import Appender 10 11from collections import defaultdict 12import warnings 13 14import numpy as np 15import pandas as pd 16from scipy import linalg as spl 17 18from statsmodels.stats.correlation_tools import cov_nearest 19from statsmodels.tools.sm_exceptions import ( 20 ConvergenceWarning, 21 NotImplementedWarning, 22 OutputWarning, 23) 24from statsmodels.tools.validation import bool_like 25 26 27class CovStruct(object): 28 """ 29 Base class for correlation and covariance structures. 30 31 An implementation of this class takes the residuals from a 32 regression model that has been fit to grouped data, and uses 33 them to estimate the within-group dependence structure of the 34 random errors in the model. 35 36 The current state of the covariance structure is represented 37 through the value of the `dep_params` attribute. 38 39 The default state of a newly-created instance should always be 40 the identity correlation matrix. 41 """ 42 43 def __init__(self, cov_nearest_method="clipped"): 44 45 # Parameters describing the dependency structure 46 self.dep_params = None 47 48 # Keep track of the number of times that the covariance was 49 # adjusted. 50 self.cov_adjust = [] 51 52 # Method for projecting the covariance matrix if it is not 53 # PSD. 54 self.cov_nearest_method = cov_nearest_method 55 56 def initialize(self, model): 57 """ 58 Called by GEE, used by implementations that need additional 59 setup prior to running `fit`. 60 61 Parameters 62 ---------- 63 model : GEE class 64 A reference to the parent GEE class instance. 65 """ 66 self.model = model 67 68 def update(self, params): 69 """ 70 Update the association parameter values based on the current 71 regression coefficients. 72 73 Parameters 74 ---------- 75 params : array_like 76 Working values for the regression parameters. 77 """ 78 raise NotImplementedError 79 80 def covariance_matrix(self, endog_expval, index): 81 """ 82 Returns the working covariance or correlation matrix for a 83 given cluster of data. 84 85 Parameters 86 ---------- 87 endog_expval : array_like 88 The expected values of endog for the cluster for which the 89 covariance or correlation matrix will be returned 90 index : int 91 The index of the cluster for which the covariance or 92 correlation matrix will be returned 93 94 Returns 95 ------- 96 M : matrix 97 The covariance or correlation matrix of endog 98 is_cor : bool 99 True if M is a correlation matrix, False if M is a 100 covariance matrix 101 """ 102 raise NotImplementedError 103 104 def covariance_matrix_solve(self, expval, index, stdev, rhs): 105 """ 106 Solves matrix equations of the form `covmat * soln = rhs` and 107 returns the values of `soln`, where `covmat` is the covariance 108 matrix represented by this class. 109 110 Parameters 111 ---------- 112 expval : array_like 113 The expected value of endog for each observed value in the 114 group. 115 index : int 116 The group index. 117 stdev : array_like 118 The standard deviation of endog for each observation in 119 the group. 120 rhs : list/tuple of array_like 121 A set of right-hand sides; each defines a matrix equation 122 to be solved. 123 124 Returns 125 ------- 126 soln : list/tuple of array_like 127 The solutions to the matrix equations. 128 129 Notes 130 ----- 131 Returns None if the solver fails. 132 133 Some dependence structures do not use `expval` and/or `index` 134 to determine the correlation matrix. Some families 135 (e.g. binomial) do not use the `stdev` parameter when forming 136 the covariance matrix. 137 138 If the covariance matrix is singular or not SPD, it is 139 projected to the nearest such matrix. These projection events 140 are recorded in the fit_history attribute of the GEE model. 141 142 Systems of linear equations with the covariance matrix as the 143 left hand side (LHS) are solved for different right hand sides 144 (RHS); the LHS is only factorized once to save time. 145 146 This is a default implementation, it can be reimplemented in 147 subclasses to optimize the linear algebra according to the 148 structure of the covariance matrix. 149 """ 150 151 vmat, is_cor = self.covariance_matrix(expval, index) 152 if is_cor: 153 vmat *= np.outer(stdev, stdev) 154 155 # Factor the covariance matrix. If the factorization fails, 156 # attempt to condition it into a factorizable matrix. 157 threshold = 1e-2 158 success = False 159 cov_adjust = 0 160 for itr in range(20): 161 try: 162 vco = spl.cho_factor(vmat) 163 success = True 164 break 165 except np.linalg.LinAlgError: 166 vmat = cov_nearest(vmat, method=self.cov_nearest_method, 167 threshold=threshold) 168 threshold *= 2 169 cov_adjust += 1 170 msg = "At least one covariance matrix was not PSD " 171 msg += "and required projection." 172 warnings.warn(msg) 173 174 self.cov_adjust.append(cov_adjust) 175 176 # Last resort if we still cannot factor the covariance matrix. 177 if not success: 178 warnings.warn( 179 "Unable to condition covariance matrix to an SPD " 180 "matrix using cov_nearest", ConvergenceWarning) 181 vmat = np.diag(np.diag(vmat)) 182 vco = spl.cho_factor(vmat) 183 184 soln = [spl.cho_solve(vco, x) for x in rhs] 185 return soln 186 187 def summary(self): 188 """ 189 Returns a text summary of the current estimate of the 190 dependence structure. 191 """ 192 raise NotImplementedError 193 194 195class Independence(CovStruct): 196 """ 197 An independence working dependence structure. 198 """ 199 200 @Appender(CovStruct.update.__doc__) 201 def update(self, params): 202 # Nothing to update 203 return 204 205 @Appender(CovStruct.covariance_matrix.__doc__) 206 def covariance_matrix(self, expval, index): 207 dim = len(expval) 208 return np.eye(dim, dtype=np.float64), True 209 210 @Appender(CovStruct.covariance_matrix_solve.__doc__) 211 def covariance_matrix_solve(self, expval, index, stdev, rhs): 212 v = stdev ** 2 213 rslt = [] 214 for x in rhs: 215 if x.ndim == 1: 216 rslt.append(x / v) 217 else: 218 rslt.append(x / v[:, None]) 219 return rslt 220 221 def summary(self): 222 return ("Observations within a cluster are modeled " 223 "as being independent.") 224 225class Unstructured(CovStruct): 226 """ 227 An unstructured dependence structure. 228 229 To use the unstructured dependence structure, a `time` 230 argument must be provided when creating the GEE. The 231 time argument must be of integer dtype, and indicates 232 which position in a complete data vector is occupied 233 by each observed value. 234 """ 235 236 def __init__(self, cov_nearest_method="clipped"): 237 238 super(Unstructured, self).__init__(cov_nearest_method) 239 240 def initialize(self, model): 241 242 self.model = model 243 244 import numbers 245 if not issubclass(self.model.time.dtype.type, numbers.Integral): 246 msg = "time must be provided and must have integer dtype" 247 raise ValueError(msg) 248 249 q = self.model.time[:, 0].max() + 1 250 251 self.dep_params = np.eye(q) 252 253 @Appender(CovStruct.covariance_matrix.__doc__) 254 def covariance_matrix(self, endog_expval, index): 255 256 if hasattr(self.model, "time"): 257 time_li = self.model.time_li 258 ix = time_li[index][:, 0] 259 return self.dep_params[np.ix_(ix, ix)],True 260 261 return self.dep_params, True 262 263 @Appender(CovStruct.update.__doc__) 264 def update(self, params): 265 266 endog = self.model.endog_li 267 nobs = self.model.nobs 268 varfunc = self.model.family.variance 269 cached_means = self.model.cached_means 270 has_weights = self.model.weights is not None 271 weights_li = self.model.weights 272 273 time_li = self.model.time_li 274 q = self.model.time.max() + 1 275 csum = np.zeros((q, q)) 276 wsum = 0. 277 cov = np.zeros((q, q)) 278 279 scale = 0. 280 for i in range(self.model.num_group): 281 282 # Get the Pearson residuals 283 expval, _ = cached_means[i] 284 stdev = np.sqrt(varfunc(expval)) 285 resid = (endog[i] - expval) / stdev 286 287 ix = time_li[i][:, 0] 288 m = np.outer(resid, resid) 289 ssr = np.sum(np.diag(m)) 290 291 w = weights_li[i] if has_weights else 1. 292 csum[np.ix_(ix, ix)] += w 293 wsum += w * len(ix) 294 cov[np.ix_(ix, ix)] += w * m 295 scale += w * ssr 296 ddof = self.model.ddof_scale 297 scale /= wsum * (nobs - ddof) / float(nobs) 298 cov /= (csum - ddof) 299 300 sd = np.sqrt(np.diag(cov)) 301 cov /= np.outer(sd, sd) 302 303 self.dep_params = cov 304 305 def summary(self): 306 print("Estimated covariance structure:") 307 print(self.dep_params) 308 309 310class Exchangeable(CovStruct): 311 """ 312 An exchangeable working dependence structure. 313 """ 314 315 def __init__(self): 316 317 super(Exchangeable, self).__init__() 318 319 # The correlation between any two values in the same cluster 320 self.dep_params = 0. 321 322 @Appender(CovStruct.update.__doc__) 323 def update(self, params): 324 325 endog = self.model.endog_li 326 327 nobs = self.model.nobs 328 329 varfunc = self.model.family.variance 330 331 cached_means = self.model.cached_means 332 333 has_weights = self.model.weights is not None 334 weights_li = self.model.weights 335 336 residsq_sum, scale = 0, 0 337 fsum1, fsum2, n_pairs = 0., 0., 0. 338 for i in range(self.model.num_group): 339 expval, _ = cached_means[i] 340 stdev = np.sqrt(varfunc(expval)) 341 resid = (endog[i] - expval) / stdev 342 f = weights_li[i] if has_weights else 1. 343 344 ssr = np.sum(resid * resid) 345 scale += f * ssr 346 fsum1 += f * len(endog[i]) 347 348 residsq_sum += f * (resid.sum() ** 2 - ssr) / 2 349 ngrp = len(resid) 350 npr = 0.5 * ngrp * (ngrp - 1) 351 fsum2 += f * npr 352 n_pairs += npr 353 354 ddof = self.model.ddof_scale 355 scale /= (fsum1 * (nobs - ddof) / float(nobs)) 356 residsq_sum /= scale 357 self.dep_params = residsq_sum / \ 358 (fsum2 * (n_pairs - ddof) / float(n_pairs)) 359 360 @Appender(CovStruct.covariance_matrix.__doc__) 361 def covariance_matrix(self, expval, index): 362 dim = len(expval) 363 dp = self.dep_params * np.ones((dim, dim), dtype=np.float64) 364 np.fill_diagonal(dp, 1) 365 return dp, True 366 367 @Appender(CovStruct.covariance_matrix_solve.__doc__) 368 def covariance_matrix_solve(self, expval, index, stdev, rhs): 369 370 k = len(expval) 371 c = self.dep_params / (1. - self.dep_params) 372 c /= 1. + self.dep_params * (k - 1) 373 374 rslt = [] 375 for x in rhs: 376 if x.ndim == 1: 377 x1 = x / stdev 378 y = x1 / (1. - self.dep_params) 379 y -= c * sum(x1) 380 y /= stdev 381 else: 382 x1 = x / stdev[:, None] 383 y = x1 / (1. - self.dep_params) 384 y -= c * x1.sum(0) 385 y /= stdev[:, None] 386 rslt.append(y) 387 388 return rslt 389 390 def summary(self): 391 return ("The correlation between two observations in the " + 392 "same cluster is %.3f" % self.dep_params) 393 394 395class Nested(CovStruct): 396 """ 397 A nested working dependence structure. 398 399 A nested working dependence structure captures unique variance 400 associated with each level in a hierarchy of partitions of the 401 cases. For each level of the hierarchy, there is a set of iid 402 random effects with mean zero, and with variance that is specific 403 to the level. These variance parameters are estimated from the 404 data using the method of moments. 405 406 The top level of the hierarchy is always defined by the required 407 `groups` argument to GEE. 408 409 The `dep_data` argument used to create the GEE defines the 410 remaining levels of the hierarchy. it should be either an array, 411 or if using the formula interface, a string that contains a 412 formula. If an array, it should contain a `n_obs x k` matrix of 413 labels, corresponding to the k levels of partitioning that are 414 nested under the top-level `groups` of the GEE instance. These 415 subgroups should be nested from left to right, so that two 416 observations with the same label for column j of `dep_data` should 417 also have the same label for all columns j' < j (this only applies 418 to observations in the same top-level cluster given by the 419 `groups` argument to GEE). 420 421 If `dep_data` is a formula, it should usually be of the form `0 + 422 a + b + ...`, where `a`, `b`, etc. contain labels defining group 423 membership. The `0 + ` should be included to prevent creation of 424 an intercept. The variable values are interpreted as labels for 425 group membership, but the variables should not be explicitly coded 426 as categorical, i.e. use `0 + a` not `0 + C(a)`. 427 428 Notes 429 ----- 430 The calculations for the nested structure involve all pairs of 431 observations within the top level `group` passed to GEE. Large 432 group sizes will result in slow iterations. 433 """ 434 435 def initialize(self, model): 436 """ 437 Called on the first call to update 438 439 `ilabels` is a list of n_i x n_i matrices containing integer 440 labels that correspond to specific correlation parameters. 441 Two elements of ilabels[i] with the same label share identical 442 variance components. 443 444 `designx` is a matrix, with each row containing dummy 445 variables indicating which variance components are associated 446 with the corresponding element of QY. 447 """ 448 449 super(Nested, self).initialize(model) 450 451 if self.model.weights is not None: 452 warnings.warn("weights not implemented for nested cov_struct, " 453 "using unweighted covariance estimate", 454 NotImplementedWarning) 455 456 # A bit of processing of the nest data 457 id_matrix = np.asarray(self.model.dep_data) 458 if id_matrix.ndim == 1: 459 id_matrix = id_matrix[:, None] 460 self.id_matrix = id_matrix 461 462 endog = self.model.endog_li 463 designx, ilabels = [], [] 464 465 # The number of layers of nesting 466 n_nest = self.id_matrix.shape[1] 467 468 for i in range(self.model.num_group): 469 ngrp = len(endog[i]) 470 glab = self.model.group_labels[i] 471 rix = self.model.group_indices[glab] 472 473 # Determine the number of common variance components 474 # shared by each pair of observations. 475 ix1, ix2 = np.tril_indices(ngrp, -1) 476 ncm = (self.id_matrix[rix[ix1], :] == 477 self.id_matrix[rix[ix2], :]).sum(1) 478 479 # This is used to construct the working correlation 480 # matrix. 481 ilabel = np.zeros((ngrp, ngrp), dtype=np.int32) 482 ilabel[(ix1, ix2)] = ncm + 1 483 ilabel[(ix2, ix1)] = ncm + 1 484 ilabels.append(ilabel) 485 486 # This is used to estimate the variance components. 487 dsx = np.zeros((len(ix1), n_nest + 1), dtype=np.float64) 488 dsx[:, 0] = 1 489 for k in np.unique(ncm): 490 ii = np.flatnonzero(ncm == k) 491 dsx[ii, 1:k + 1] = 1 492 designx.append(dsx) 493 494 self.designx = np.concatenate(designx, axis=0) 495 self.ilabels = ilabels 496 497 svd = np.linalg.svd(self.designx, 0) 498 self.designx_u = svd[0] 499 self.designx_s = svd[1] 500 self.designx_v = svd[2].T 501 502 @Appender(CovStruct.update.__doc__) 503 def update(self, params): 504 505 endog = self.model.endog_li 506 507 nobs = self.model.nobs 508 dim = len(params) 509 510 if self.designx is None: 511 self._compute_design(self.model) 512 513 cached_means = self.model.cached_means 514 515 varfunc = self.model.family.variance 516 517 dvmat = [] 518 scale = 0. 519 for i in range(self.model.num_group): 520 521 expval, _ = cached_means[i] 522 523 stdev = np.sqrt(varfunc(expval)) 524 resid = (endog[i] - expval) / stdev 525 526 ix1, ix2 = np.tril_indices(len(resid), -1) 527 dvmat.append(resid[ix1] * resid[ix2]) 528 529 scale += np.sum(resid ** 2) 530 531 dvmat = np.concatenate(dvmat) 532 scale /= (nobs - dim) 533 534 # Use least squares regression to estimate the variance 535 # components 536 vcomp_coeff = np.dot(self.designx_v, np.dot(self.designx_u.T, 537 dvmat) / self.designx_s) 538 539 self.vcomp_coeff = np.clip(vcomp_coeff, 0, np.inf) 540 self.scale = scale 541 542 self.dep_params = self.vcomp_coeff.copy() 543 544 @Appender(CovStruct.covariance_matrix.__doc__) 545 def covariance_matrix(self, expval, index): 546 547 dim = len(expval) 548 549 # First iteration 550 if self.dep_params is None: 551 return np.eye(dim, dtype=np.float64), True 552 553 ilabel = self.ilabels[index] 554 555 c = np.r_[self.scale, np.cumsum(self.vcomp_coeff)] 556 vmat = c[ilabel] 557 vmat /= self.scale 558 return vmat, True 559 560 def summary(self): 561 """ 562 Returns a summary string describing the state of the 563 dependence structure. 564 """ 565 566 dep_names = ["Groups"] 567 if hasattr(self.model, "_dep_data_names"): 568 dep_names.extend(self.model._dep_data_names) 569 else: 570 dep_names.extend(["Component %d:" % (k + 1) for k in range(len(self.vcomp_coeff) - 1)]) 571 if hasattr(self.model, "_groups_name"): 572 dep_names[0] = self.model._groups_name 573 dep_names.append("Residual") 574 575 vc = self.vcomp_coeff.tolist() 576 vc.append(self.scale - np.sum(vc)) 577 578 smry = pd.DataFrame({"Variance": vc}, index=dep_names) 579 580 return smry 581 582 583class Stationary(CovStruct): 584 """ 585 A stationary covariance structure. 586 587 The correlation between two observations is an arbitrary function 588 of the distance between them. Distances up to a given maximum 589 value are included in the covariance model. 590 591 Parameters 592 ---------- 593 max_lag : float 594 The largest distance that is included in the covariance model. 595 grid : bool 596 If True, the index positions in the data (after dropping missing 597 values) are used to define distances, and the `time` variable is 598 ignored. 599 """ 600 601 def __init__(self, max_lag=1, grid=None): 602 603 super(Stationary, self).__init__() 604 grid = bool_like(grid, "grid", optional=True) 605 if grid is None: 606 warnings.warn( 607 "grid=True will become default in a future version", 608 FutureWarning 609 ) 610 611 self.max_lag = max_lag 612 self.grid = bool(grid) 613 self.dep_params = np.zeros(max_lag + 1) 614 615 def initialize(self, model): 616 617 super(Stationary, self).initialize(model) 618 619 # Time used as an index needs to be integer type. 620 if not self.grid: 621 time = self.model.time[:, 0].astype(np.int32) 622 self.time = self.model.cluster_list(time) 623 624 @Appender(CovStruct.update.__doc__) 625 def update(self, params): 626 627 if self.grid: 628 self.update_grid(params) 629 else: 630 self.update_nogrid(params) 631 632 def update_grid(self, params): 633 634 endog = self.model.endog_li 635 cached_means = self.model.cached_means 636 varfunc = self.model.family.variance 637 638 dep_params = np.zeros(self.max_lag + 1) 639 for i in range(self.model.num_group): 640 641 expval, _ = cached_means[i] 642 stdev = np.sqrt(varfunc(expval)) 643 resid = (endog[i] - expval) / stdev 644 645 dep_params[0] += np.sum(resid * resid) / len(resid) 646 for j in range(1, self.max_lag + 1): 647 v = resid[j:] 648 dep_params[j] += np.sum(resid[0:-j] * v) / len(v) 649 650 dep_params /= dep_params[0] 651 self.dep_params = dep_params 652 653 def update_nogrid(self, params): 654 655 endog = self.model.endog_li 656 cached_means = self.model.cached_means 657 varfunc = self.model.family.variance 658 659 dep_params = np.zeros(self.max_lag + 1) 660 dn = np.zeros(self.max_lag + 1) 661 resid_ssq = 0 662 resid_ssq_n = 0 663 for i in range(self.model.num_group): 664 665 expval, _ = cached_means[i] 666 stdev = np.sqrt(varfunc(expval)) 667 resid = (endog[i] - expval) / stdev 668 669 j1, j2 = np.tril_indices(len(expval), -1) 670 dx = np.abs(self.time[i][j1] - self.time[i][j2]) 671 ii = np.flatnonzero(dx <= self.max_lag) 672 j1 = j1[ii] 673 j2 = j2[ii] 674 dx = dx[ii] 675 676 vs = np.bincount(dx, weights=resid[j1] * resid[j2], 677 minlength=self.max_lag + 1) 678 vd = np.bincount(dx, minlength=self.max_lag + 1) 679 680 resid_ssq += np.sum(resid**2) 681 resid_ssq_n += len(resid) 682 683 ii = np.flatnonzero(vd > 0) 684 if len(ii) > 0: 685 dn[ii] += 1 686 dep_params[ii] += vs[ii] / vd[ii] 687 688 i0 = np.flatnonzero(dn > 0) 689 dep_params[i0] /= dn[i0] 690 resid_msq = resid_ssq / resid_ssq_n 691 dep_params /= resid_msq 692 self.dep_params = dep_params 693 694 @Appender(CovStruct.covariance_matrix.__doc__) 695 def covariance_matrix(self, endog_expval, index): 696 697 if self.grid: 698 return self.covariance_matrix_grid(endog_expval, index) 699 700 j1, j2 = np.tril_indices(len(endog_expval), -1) 701 dx = np.abs(self.time[index][j1] - self.time[index][j2]) 702 ii = np.flatnonzero(dx <= self.max_lag) 703 j1 = j1[ii] 704 j2 = j2[ii] 705 dx = dx[ii] 706 707 cmat = np.eye(len(endog_expval)) 708 cmat[j1, j2] = self.dep_params[dx] 709 cmat[j2, j1] = self.dep_params[dx] 710 711 return cmat, True 712 713 def covariance_matrix_grid(self, endog_expval, index): 714 715 from scipy.linalg import toeplitz 716 r = np.zeros(len(endog_expval)) 717 r[0] = 1 718 r[1:self.max_lag + 1] = self.dep_params[1:] 719 return toeplitz(r), True 720 721 @Appender(CovStruct.covariance_matrix_solve.__doc__) 722 def covariance_matrix_solve(self, expval, index, stdev, rhs): 723 724 if not self.grid: 725 return super(Stationary, self).covariance_matrix_solve( 726 expval, index, stdev, rhs) 727 728 from statsmodels.tools.linalg import stationary_solve 729 r = np.zeros(len(expval)) 730 r[0:self.max_lag] = self.dep_params[1:] 731 732 rslt = [] 733 for x in rhs: 734 if x.ndim == 1: 735 y = x / stdev 736 rslt.append(stationary_solve(r, y) / stdev) 737 else: 738 y = x / stdev[:, None] 739 rslt.append(stationary_solve(r, y) / stdev[:, None]) 740 741 return rslt 742 743 def summary(self): 744 745 lag = np.arange(self.max_lag + 1) 746 return pd.DataFrame({"Lag": lag, "Cov": self.dep_params}) 747 748 749class Autoregressive(CovStruct): 750 """ 751 A first-order autoregressive working dependence structure. 752 753 The dependence is defined in terms of the `time` component of the 754 parent GEE class, which defaults to the index position of each 755 value within its cluster, based on the order of values in the 756 input data set. Time represents a potentially multidimensional 757 index from which distances between pairs of observations can be 758 determined. 759 760 The correlation between two observations in the same cluster is 761 dep_params^distance, where `dep_params` contains the (scalar) 762 autocorrelation parameter to be estimated, and `distance` is the 763 distance between the two observations, calculated from their 764 corresponding time values. `time` is stored as an n_obs x k 765 matrix, where `k` represents the number of dimensions in the time 766 index. 767 768 The autocorrelation parameter is estimated using weighted 769 nonlinear least squares, regressing each value within a cluster on 770 each preceding value in the same cluster. 771 772 Parameters 773 ---------- 774 dist_func : function from R^k x R^k to R^+, optional 775 A function that computes the distance between the two 776 observations based on their `time` values. 777 778 References 779 ---------- 780 B Rosner, A Munoz. Autoregressive modeling for the analysis of 781 longitudinal data with unequally spaced examinations. Statistics 782 in medicine. Vol 7, 59-71, 1988. 783 """ 784 785 def __init__(self, dist_func=None, grid=None): 786 787 super(Autoregressive, self).__init__() 788 grid = bool_like(grid, "grid", optional=True) 789 # The function for determining distances based on time 790 if dist_func is None: 791 self.dist_func = lambda x, y: np.abs(x - y).sum() 792 else: 793 self.dist_func = dist_func 794 795 if grid is None: 796 warnings.warn( 797 "grid=True will become default in a future version", 798 FutureWarning 799 ) 800 self.grid = bool(grid) 801 if not self.grid: 802 self.designx = None 803 804 # The autocorrelation parameter 805 self.dep_params = 0. 806 807 @Appender(CovStruct.update.__doc__) 808 def update(self, params): 809 810 if self.model.weights is not None: 811 warnings.warn("weights not implemented for autoregressive " 812 "cov_struct, using unweighted covariance estimate", 813 NotImplementedWarning) 814 815 if self.grid: 816 self._update_grid(params) 817 else: 818 self._update_nogrid(params) 819 820 def _update_grid(self, params): 821 822 cached_means = self.model.cached_means 823 scale = self.model.estimate_scale() 824 varfunc = self.model.family.variance 825 endog = self.model.endog_li 826 827 lag0, lag1 = 0.0, 0.0 828 for i in range(self.model.num_group): 829 830 expval, _ = cached_means[i] 831 stdev = np.sqrt(scale * varfunc(expval)) 832 resid = (endog[i] - expval) / stdev 833 834 n = len(resid) 835 if n > 1: 836 lag1 += np.sum(resid[0:-1] * resid[1:]) / (n - 1) 837 lag0 += np.sum(resid**2) / n 838 839 self.dep_params = lag1 / lag0 840 841 def _update_nogrid(self, params): 842 843 endog = self.model.endog_li 844 time = self.model.time_li 845 846 # Only need to compute this once 847 if self.designx is not None: 848 designx = self.designx 849 else: 850 designx = [] 851 for i in range(self.model.num_group): 852 853 ngrp = len(endog[i]) 854 if ngrp == 0: 855 continue 856 857 # Loop over pairs of observations within a cluster 858 for j1 in range(ngrp): 859 for j2 in range(j1): 860 designx.append(self.dist_func(time[i][j1, :], 861 time[i][j2, :])) 862 863 designx = np.array(designx) 864 self.designx = designx 865 866 scale = self.model.estimate_scale() 867 varfunc = self.model.family.variance 868 cached_means = self.model.cached_means 869 870 # Weights 871 var = 1. - self.dep_params ** (2 * designx) 872 var /= 1. - self.dep_params ** 2 873 wts = 1. / var 874 wts /= wts.sum() 875 876 residmat = [] 877 for i in range(self.model.num_group): 878 879 expval, _ = cached_means[i] 880 stdev = np.sqrt(scale * varfunc(expval)) 881 resid = (endog[i] - expval) / stdev 882 883 ngrp = len(resid) 884 for j1 in range(ngrp): 885 for j2 in range(j1): 886 residmat.append([resid[j1], resid[j2]]) 887 888 residmat = np.array(residmat) 889 890 # Need to minimize this 891 def fitfunc(a): 892 dif = residmat[:, 0] - (a ** designx) * residmat[:, 1] 893 return np.dot(dif ** 2, wts) 894 895 # Left bracket point 896 b_lft, f_lft = 0., fitfunc(0.) 897 898 # Center bracket point 899 b_ctr, f_ctr = 0.5, fitfunc(0.5) 900 while f_ctr > f_lft: 901 b_ctr /= 2 902 f_ctr = fitfunc(b_ctr) 903 if b_ctr < 1e-8: 904 self.dep_params = 0 905 return 906 907 # Right bracket point 908 b_rgt, f_rgt = 0.75, fitfunc(0.75) 909 while f_rgt < f_ctr: 910 b_rgt = b_rgt + (1. - b_rgt) / 2 911 f_rgt = fitfunc(b_rgt) 912 if b_rgt > 1. - 1e-6: 913 raise ValueError( 914 "Autoregressive: unable to find right bracket") 915 916 from scipy.optimize import brent 917 self.dep_params = brent(fitfunc, brack=[b_lft, b_ctr, b_rgt]) 918 919 @Appender(CovStruct.covariance_matrix.__doc__) 920 def covariance_matrix(self, endog_expval, index): 921 ngrp = len(endog_expval) 922 if self.dep_params == 0: 923 return np.eye(ngrp, dtype=np.float64), True 924 idx = np.arange(ngrp) 925 cmat = self.dep_params ** np.abs(idx[:, None] - idx[None, :]) 926 return cmat, True 927 928 @Appender(CovStruct.covariance_matrix_solve.__doc__) 929 def covariance_matrix_solve(self, expval, index, stdev, rhs): 930 # The inverse of an AR(1) covariance matrix is tri-diagonal. 931 932 k = len(expval) 933 r = self.dep_params 934 soln = [] 935 936 # RHS has 1 row 937 if k == 1: 938 return [x / stdev ** 2 for x in rhs] 939 940 # RHS has 2 rows 941 if k == 2: 942 mat = np.array([[1, -r], [-r, 1]]) 943 mat /= (1. - r ** 2) 944 for x in rhs: 945 if x.ndim == 1: 946 x1 = x / stdev 947 else: 948 x1 = x / stdev[:, None] 949 x1 = np.dot(mat, x1) 950 if x.ndim == 1: 951 x1 /= stdev 952 else: 953 x1 /= stdev[:, None] 954 soln.append(x1) 955 return soln 956 957 # RHS has >= 3 rows: values c0, c1, c2 defined below give 958 # the inverse. c0 is on the diagonal, except for the first 959 # and last position. c1 is on the first and last position of 960 # the diagonal. c2 is on the sub/super diagonal. 961 c0 = (1. + r ** 2) / (1. - r ** 2) 962 c1 = 1. / (1. - r ** 2) 963 c2 = -r / (1. - r ** 2) 964 soln = [] 965 for x in rhs: 966 flatten = False 967 if x.ndim == 1: 968 x = x[:, None] 969 flatten = True 970 x1 = x / stdev[:, None] 971 972 z0 = np.zeros((1, x1.shape[1])) 973 rhs1 = np.concatenate((x1[1:, :], z0), axis=0) 974 rhs2 = np.concatenate((z0, x1[0:-1, :]), axis=0) 975 976 y = c0 * x1 + c2 * rhs1 + c2 * rhs2 977 y[0, :] = c1 * x1[0, :] + c2 * x1[1, :] 978 y[-1, :] = c1 * x1[-1, :] + c2 * x1[-2, :] 979 980 y /= stdev[:, None] 981 982 if flatten: 983 y = np.squeeze(y) 984 985 soln.append(y) 986 987 return soln 988 989 def summary(self): 990 991 return ("Autoregressive(1) dependence parameter: %.3f\n" % 992 self.dep_params) 993 994 995class CategoricalCovStruct(CovStruct): 996 """ 997 Parent class for covariance structure for categorical data models. 998 999 Attributes 1000 ---------- 1001 nlevel : int 1002 The number of distinct levels for the outcome variable. 1003 ibd : list 1004 A list whose i^th element ibd[i] is an array whose rows 1005 contain integer pairs (a,b), where endog_li[i][a:b] is the 1006 subvector of binary indicators derived from the same ordinal 1007 value. 1008 """ 1009 1010 def initialize(self, model): 1011 1012 super(CategoricalCovStruct, self).initialize(model) 1013 1014 self.nlevel = len(model.endog_values) 1015 self._ncut = self.nlevel - 1 1016 1017 from numpy.lib.stride_tricks import as_strided 1018 b = np.dtype(np.int64).itemsize 1019 1020 ibd = [] 1021 for v in model.endog_li: 1022 jj = np.arange(0, len(v) + 1, self._ncut, dtype=np.int64) 1023 jj = as_strided(jj, shape=(len(jj) - 1, 2), strides=(b, b)) 1024 ibd.append(jj) 1025 1026 self.ibd = ibd 1027 1028 1029class GlobalOddsRatio(CategoricalCovStruct): 1030 """ 1031 Estimate the global odds ratio for a GEE with ordinal or nominal 1032 data. 1033 1034 References 1035 ---------- 1036 PJ Heagerty and S Zeger. "Marginal Regression Models for Clustered 1037 Ordinal Measurements". Journal of the American Statistical 1038 Association Vol. 91, Issue 435 (1996). 1039 1040 Thomas Lumley. Generalized Estimating Equations for Ordinal Data: 1041 A Note on Working Correlation Structures. Biometrics Vol. 52, 1042 No. 1 (Mar., 1996), pp. 354-361 1043 http://www.jstor.org/stable/2533173 1044 1045 Notes 1046 ----- 1047 The following data structures are calculated in the class: 1048 1049 'ibd' is a list whose i^th element ibd[i] is a sequence of integer 1050 pairs (a,b), where endog_li[i][a:b] is the subvector of binary 1051 indicators derived from the same ordinal value. 1052 1053 `cpp` is a dictionary where cpp[group] is a map from cut-point 1054 pairs (c,c') to the indices of all between-subject pairs derived 1055 from the given cut points. 1056 """ 1057 1058 def __init__(self, endog_type): 1059 super(GlobalOddsRatio, self).__init__() 1060 self.endog_type = endog_type 1061 self.dep_params = 0. 1062 1063 def initialize(self, model): 1064 1065 super(GlobalOddsRatio, self).initialize(model) 1066 1067 if self.model.weights is not None: 1068 warnings.warn("weights not implemented for GlobalOddsRatio " 1069 "cov_struct, using unweighted covariance estimate", 1070 NotImplementedWarning) 1071 1072 # Need to restrict to between-subject pairs 1073 cpp = [] 1074 for v in model.endog_li: 1075 1076 # Number of subjects in this group 1077 m = int(len(v) / self._ncut) 1078 i1, i2 = np.tril_indices(m, -1) 1079 1080 cpp1 = {} 1081 for k1 in range(self._ncut): 1082 for k2 in range(k1 + 1): 1083 jj = np.zeros((len(i1), 2), dtype=np.int64) 1084 jj[:, 0] = i1 * self._ncut + k1 1085 jj[:, 1] = i2 * self._ncut + k2 1086 cpp1[(k2, k1)] = jj 1087 1088 cpp.append(cpp1) 1089 1090 self.cpp = cpp 1091 1092 # Initialize the dependence parameters 1093 self.crude_or = self.observed_crude_oddsratio() 1094 if self.model.update_dep: 1095 self.dep_params = self.crude_or 1096 1097 def pooled_odds_ratio(self, tables): 1098 """ 1099 Returns the pooled odds ratio for a list of 2x2 tables. 1100 1101 The pooled odds ratio is the inverse variance weighted average 1102 of the sample odds ratios of the tables. 1103 """ 1104 1105 if len(tables) == 0: 1106 return 1. 1107 1108 # Get the sampled odds ratios and variances 1109 log_oddsratio, var = [], [] 1110 for table in tables: 1111 lor = np.log(table[1, 1]) + np.log(table[0, 0]) -\ 1112 np.log(table[0, 1]) - np.log(table[1, 0]) 1113 log_oddsratio.append(lor) 1114 var.append((1 / table.astype(np.float64)).sum()) 1115 1116 # Calculate the inverse variance weighted average 1117 wts = [1 / v for v in var] 1118 wtsum = sum(wts) 1119 wts = [w / wtsum for w in wts] 1120 log_pooled_or = sum([w * e for w, e in zip(wts, log_oddsratio)]) 1121 1122 return np.exp(log_pooled_or) 1123 1124 @Appender(CovStruct.covariance_matrix.__doc__) 1125 def covariance_matrix(self, expected_value, index): 1126 1127 vmat = self.get_eyy(expected_value, index) 1128 vmat -= np.outer(expected_value, expected_value) 1129 return vmat, False 1130 1131 def observed_crude_oddsratio(self): 1132 """ 1133 To obtain the crude (global) odds ratio, first pool all binary 1134 indicators corresponding to a given pair of cut points (c,c'), 1135 then calculate the odds ratio for this 2x2 table. The crude 1136 odds ratio is the inverse variance weighted average of these 1137 odds ratios. Since the covariate effects are ignored, this OR 1138 will generally be greater than the stratified OR. 1139 """ 1140 1141 cpp = self.cpp 1142 endog = self.model.endog_li 1143 1144 # Storage for the contingency tables for each (c,c') 1145 tables = {} 1146 for ii in cpp[0].keys(): 1147 tables[ii] = np.zeros((2, 2), dtype=np.float64) 1148 1149 # Get the observed crude OR 1150 for i in range(len(endog)): 1151 1152 # The observed joint values for the current cluster 1153 yvec = endog[i] 1154 endog_11 = np.outer(yvec, yvec) 1155 endog_10 = np.outer(yvec, 1. - yvec) 1156 endog_01 = np.outer(1. - yvec, yvec) 1157 endog_00 = np.outer(1. - yvec, 1. - yvec) 1158 1159 cpp1 = cpp[i] 1160 for ky in cpp1.keys(): 1161 ix = cpp1[ky] 1162 tables[ky][1, 1] += endog_11[ix[:, 0], ix[:, 1]].sum() 1163 tables[ky][1, 0] += endog_10[ix[:, 0], ix[:, 1]].sum() 1164 tables[ky][0, 1] += endog_01[ix[:, 0], ix[:, 1]].sum() 1165 tables[ky][0, 0] += endog_00[ix[:, 0], ix[:, 1]].sum() 1166 1167 return self.pooled_odds_ratio(list(tables.values())) 1168 1169 def get_eyy(self, endog_expval, index): 1170 """ 1171 Returns a matrix V such that V[i,j] is the joint probability 1172 that endog[i] = 1 and endog[j] = 1, based on the marginal 1173 probabilities of endog and the global odds ratio `current_or`. 1174 """ 1175 1176 current_or = self.dep_params 1177 ibd = self.ibd[index] 1178 1179 # The between-observation joint probabilities 1180 if current_or == 1.0: 1181 vmat = np.outer(endog_expval, endog_expval) 1182 else: 1183 psum = endog_expval[:, None] + endog_expval[None, :] 1184 pprod = endog_expval[:, None] * endog_expval[None, :] 1185 pfac = np.sqrt((1. + psum * (current_or - 1.)) ** 2 + 1186 4 * current_or * (1. - current_or) * pprod) 1187 vmat = 1. + psum * (current_or - 1.) - pfac 1188 vmat /= 2. * (current_or - 1) 1189 1190 # Fix E[YY'] for elements that belong to same observation 1191 for bdl in ibd: 1192 evy = endog_expval[bdl[0]:bdl[1]] 1193 if self.endog_type == "ordinal": 1194 vmat[bdl[0]:bdl[1], bdl[0]:bdl[1]] =\ 1195 np.minimum.outer(evy, evy) 1196 else: 1197 vmat[bdl[0]:bdl[1], bdl[0]:bdl[1]] = np.diag(evy) 1198 1199 return vmat 1200 1201 @Appender(CovStruct.update.__doc__) 1202 def update(self, params): 1203 """ 1204 Update the global odds ratio based on the current value of 1205 params. 1206 """ 1207 1208 cpp = self.cpp 1209 cached_means = self.model.cached_means 1210 1211 # This will happen if all the clusters have only 1212 # one observation 1213 if len(cpp[0]) == 0: 1214 return 1215 1216 tables = {} 1217 for ii in cpp[0]: 1218 tables[ii] = np.zeros((2, 2), dtype=np.float64) 1219 1220 for i in range(self.model.num_group): 1221 1222 endog_expval, _ = cached_means[i] 1223 1224 emat_11 = self.get_eyy(endog_expval, i) 1225 emat_10 = endog_expval[:, None] - emat_11 1226 emat_01 = -emat_11 + endog_expval 1227 emat_00 = 1. - (emat_11 + emat_10 + emat_01) 1228 1229 cpp1 = cpp[i] 1230 for ky in cpp1.keys(): 1231 ix = cpp1[ky] 1232 tables[ky][1, 1] += emat_11[ix[:, 0], ix[:, 1]].sum() 1233 tables[ky][1, 0] += emat_10[ix[:, 0], ix[:, 1]].sum() 1234 tables[ky][0, 1] += emat_01[ix[:, 0], ix[:, 1]].sum() 1235 tables[ky][0, 0] += emat_00[ix[:, 0], ix[:, 1]].sum() 1236 1237 cor_expval = self.pooled_odds_ratio(list(tables.values())) 1238 1239 self.dep_params *= self.crude_or / cor_expval 1240 if not np.isfinite(self.dep_params): 1241 self.dep_params = 1. 1242 warnings.warn("dep_params became inf, resetting to 1", 1243 ConvergenceWarning) 1244 1245 def summary(self): 1246 return "Global odds ratio: %.3f\n" % self.dep_params 1247 1248 1249class OrdinalIndependence(CategoricalCovStruct): 1250 """ 1251 An independence covariance structure for ordinal models. 1252 1253 The working covariance between indicators derived from different 1254 observations is zero. The working covariance between indicators 1255 derived form a common observation is determined from their current 1256 mean values. 1257 1258 There are no parameters to estimate in this covariance structure. 1259 """ 1260 1261 def covariance_matrix(self, expected_value, index): 1262 1263 ibd = self.ibd[index] 1264 n = len(expected_value) 1265 vmat = np.zeros((n, n)) 1266 1267 for bdl in ibd: 1268 ev = expected_value[bdl[0]:bdl[1]] 1269 vmat[bdl[0]:bdl[1], bdl[0]:bdl[1]] =\ 1270 np.minimum.outer(ev, ev) - np.outer(ev, ev) 1271 1272 return vmat, False 1273 1274 # Nothing to update 1275 def update(self, params): 1276 pass 1277 1278 1279class NominalIndependence(CategoricalCovStruct): 1280 """ 1281 An independence covariance structure for nominal models. 1282 1283 The working covariance between indicators derived from different 1284 observations is zero. The working covariance between indicators 1285 derived form a common observation is determined from their current 1286 mean values. 1287 1288 There are no parameters to estimate in this covariance structure. 1289 """ 1290 1291 def covariance_matrix(self, expected_value, index): 1292 1293 ibd = self.ibd[index] 1294 n = len(expected_value) 1295 vmat = np.zeros((n, n)) 1296 1297 for bdl in ibd: 1298 ev = expected_value[bdl[0]:bdl[1]] 1299 vmat[bdl[0]:bdl[1], bdl[0]:bdl[1]] =\ 1300 np.diag(ev) - np.outer(ev, ev) 1301 1302 return vmat, False 1303 1304 # Nothing to update 1305 def update(self, params): 1306 pass 1307 1308 1309class Equivalence(CovStruct): 1310 """ 1311 A covariance structure defined in terms of equivalence classes. 1312 1313 An 'equivalence class' is a set of pairs of observations such that 1314 the covariance of every pair within the equivalence class has a 1315 common value. 1316 1317 Parameters 1318 ---------- 1319 pairs : dict-like 1320 A dictionary of dictionaries, where `pairs[group][label]` 1321 provides the indices of all pairs of observations in the group 1322 that have the same covariance value. Specifically, 1323 `pairs[group][label]` is a tuple `(j1, j2)`, where `j1` and `j2` 1324 are integer arrays of the same length. `j1[i], j2[i]` is one 1325 index pair that belongs to the `label` equivalence class. Only 1326 one triangle of each covariance matrix should be included. 1327 Positions where j1 and j2 have the same value are variance 1328 parameters. 1329 labels : array_like 1330 An array of labels such that every distinct pair of labels 1331 defines an equivalence class. Either `labels` or `pairs` must 1332 be provided. When the two labels in a pair are equal two 1333 equivalence classes are defined: one for the diagonal elements 1334 (corresponding to variances) and one for the off-diagonal 1335 elements (corresponding to covariances). 1336 return_cov : bool 1337 If True, `covariance_matrix` returns an estimate of the 1338 covariance matrix, otherwise returns an estimate of the 1339 correlation matrix. 1340 1341 Notes 1342 ----- 1343 Using `labels` to define the class is much easier than using 1344 `pairs`, but is less general. 1345 1346 Any pair of values not contained in `pairs` will be assigned zero 1347 covariance. 1348 1349 The index values in `pairs` are row indices into the `exog` 1350 matrix. They are not updated if missing data are present. When 1351 using this covariance structure, missing data should be removed 1352 before constructing the model. 1353 1354 If using `labels`, after a model is defined using the covariance 1355 structure it is possible to remove a label pair from the second 1356 level of the `pairs` dictionary to force the corresponding 1357 covariance to be zero. 1358 1359 Examples 1360 -------- 1361 The following sets up the `pairs` dictionary for a model with two 1362 groups, equal variance for all observations, and constant 1363 covariance for all pairs of observations within each group. 1364 1365 >> pairs = {0: {}, 1: {}} 1366 >> pairs[0][0] = (np.r_[0, 1, 2], np.r_[0, 1, 2]) 1367 >> pairs[0][1] = np.tril_indices(3, -1) 1368 >> pairs[1][0] = (np.r_[3, 4, 5], np.r_[3, 4, 5]) 1369 >> pairs[1][2] = 3 + np.tril_indices(3, -1) 1370 """ 1371 1372 def __init__(self, pairs=None, labels=None, return_cov=False): 1373 1374 super(Equivalence, self).__init__() 1375 1376 if (pairs is None) and (labels is None): 1377 raise ValueError( 1378 "Equivalence cov_struct requires either `pairs` or `labels`") 1379 1380 if (pairs is not None) and (labels is not None): 1381 raise ValueError( 1382 "Equivalence cov_struct accepts only one of `pairs` " 1383 "and `labels`") 1384 1385 if pairs is not None: 1386 import copy 1387 self.pairs = copy.deepcopy(pairs) 1388 1389 if labels is not None: 1390 self.labels = np.asarray(labels) 1391 1392 self.return_cov = return_cov 1393 1394 def _make_pairs(self, i, j): 1395 """ 1396 Create arrays containing all unique ordered pairs of i, j. 1397 1398 The arrays i and j must be one-dimensional containing non-negative 1399 integers. 1400 """ 1401 1402 mat = np.zeros((len(i) * len(j), 2), dtype=np.int32) 1403 1404 # Create the pairs and order them 1405 f = np.ones(len(j)) 1406 mat[:, 0] = np.kron(f, i).astype(np.int32) 1407 f = np.ones(len(i)) 1408 mat[:, 1] = np.kron(j, f).astype(np.int32) 1409 mat.sort(1) 1410 1411 # Remove repeated rows 1412 try: 1413 dtype = np.dtype((np.void, mat.dtype.itemsize * mat.shape[1])) 1414 bmat = np.ascontiguousarray(mat).view(dtype) 1415 _, idx = np.unique(bmat, return_index=True) 1416 except TypeError: 1417 # workaround for old numpy that cannot call unique with complex 1418 # dtypes 1419 rs = np.random.RandomState(4234) 1420 bmat = np.dot(mat, rs.uniform(size=mat.shape[1])) 1421 _, idx = np.unique(bmat, return_index=True) 1422 mat = mat[idx, :] 1423 1424 return mat[:, 0], mat[:, 1] 1425 1426 def _pairs_from_labels(self): 1427 1428 from collections import defaultdict 1429 pairs = defaultdict(lambda: defaultdict(lambda: None)) 1430 1431 model = self.model 1432 1433 df = pd.DataFrame({"labels": self.labels, "groups": model.groups}) 1434 gb = df.groupby(["groups", "labels"]) 1435 1436 ulabels = np.unique(self.labels) 1437 1438 for g_ix, g_lb in enumerate(model.group_labels): 1439 1440 # Loop over label pairs 1441 for lx1 in range(len(ulabels)): 1442 for lx2 in range(lx1 + 1): 1443 1444 lb1 = ulabels[lx1] 1445 lb2 = ulabels[lx2] 1446 1447 try: 1448 i1 = gb.groups[(g_lb, lb1)] 1449 i2 = gb.groups[(g_lb, lb2)] 1450 except KeyError: 1451 continue 1452 1453 i1, i2 = self._make_pairs(i1, i2) 1454 1455 clabel = str(lb1) + "/" + str(lb2) 1456 1457 # Variance parameters belong in their own equiv class. 1458 jj = np.flatnonzero(i1 == i2) 1459 if len(jj) > 0: 1460 clabelv = clabel + "/v" 1461 pairs[g_lb][clabelv] = (i1[jj], i2[jj]) 1462 1463 # Covariance parameters 1464 jj = np.flatnonzero(i1 != i2) 1465 if len(jj) > 0: 1466 i1 = i1[jj] 1467 i2 = i2[jj] 1468 pairs[g_lb][clabel] = (i1, i2) 1469 1470 self.pairs = pairs 1471 1472 def initialize(self, model): 1473 1474 super(Equivalence, self).initialize(model) 1475 1476 if self.model.weights is not None: 1477 warnings.warn("weights not implemented for equalence cov_struct, " 1478 "using unweighted covariance estimate", 1479 NotImplementedWarning) 1480 1481 if not hasattr(self, 'pairs'): 1482 self._pairs_from_labels() 1483 1484 # Initialize so that any equivalence class containing a 1485 # variance parameter has value 1. 1486 self.dep_params = defaultdict(lambda: 0.) 1487 self._var_classes = set([]) 1488 for gp in self.model.group_labels: 1489 for lb in self.pairs[gp]: 1490 j1, j2 = self.pairs[gp][lb] 1491 if np.any(j1 == j2): 1492 if not np.all(j1 == j2): 1493 warnings.warn( 1494 "equivalence class contains both variance " 1495 "and covariance parameters", OutputWarning) 1496 self._var_classes.add(lb) 1497 self.dep_params[lb] = 1 1498 1499 # Need to start indexing at 0 within each group. 1500 # rx maps olds indices to new indices 1501 rx = -1 * np.ones(len(self.model.endog), dtype=np.int32) 1502 for g_ix, g_lb in enumerate(self.model.group_labels): 1503 ii = self.model.group_indices[g_lb] 1504 rx[ii] = np.arange(len(ii), dtype=np.int32) 1505 1506 # Reindex 1507 for gp in self.model.group_labels: 1508 for lb in self.pairs[gp].keys(): 1509 a, b = self.pairs[gp][lb] 1510 self.pairs[gp][lb] = (rx[a], rx[b]) 1511 1512 @Appender(CovStruct.update.__doc__) 1513 def update(self, params): 1514 1515 endog = self.model.endog_li 1516 varfunc = self.model.family.variance 1517 cached_means = self.model.cached_means 1518 dep_params = defaultdict(lambda: [0., 0., 0.]) 1519 n_pairs = defaultdict(lambda: 0) 1520 dim = len(params) 1521 1522 for k, gp in enumerate(self.model.group_labels): 1523 expval, _ = cached_means[k] 1524 stdev = np.sqrt(varfunc(expval)) 1525 resid = (endog[k] - expval) / stdev 1526 for lb in self.pairs[gp].keys(): 1527 if (not self.return_cov) and lb in self._var_classes: 1528 continue 1529 jj = self.pairs[gp][lb] 1530 dep_params[lb][0] += np.sum(resid[jj[0]] * resid[jj[1]]) 1531 if not self.return_cov: 1532 dep_params[lb][1] += np.sum(resid[jj[0]] ** 2) 1533 dep_params[lb][2] += np.sum(resid[jj[1]] ** 2) 1534 n_pairs[lb] += len(jj[0]) 1535 1536 if self.return_cov: 1537 for lb in dep_params.keys(): 1538 dep_params[lb] = dep_params[lb][0] / (n_pairs[lb] - dim) 1539 else: 1540 for lb in dep_params.keys(): 1541 den = np.sqrt(dep_params[lb][1] * dep_params[lb][2]) 1542 dep_params[lb] = dep_params[lb][0] / den 1543 for lb in self._var_classes: 1544 dep_params[lb] = 1. 1545 1546 self.dep_params = dep_params 1547 self.n_pairs = n_pairs 1548 1549 @Appender(CovStruct.covariance_matrix.__doc__) 1550 def covariance_matrix(self, expval, index): 1551 dim = len(expval) 1552 cmat = np.zeros((dim, dim)) 1553 g_lb = self.model.group_labels[index] 1554 1555 for lb in self.pairs[g_lb].keys(): 1556 j1, j2 = self.pairs[g_lb][lb] 1557 cmat[j1, j2] = self.dep_params[lb] 1558 1559 cmat = cmat + cmat.T 1560 np.fill_diagonal(cmat, cmat.diagonal() / 2) 1561 1562 return cmat, not self.return_cov 1563