Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1from statsmodels.compat.python import lzip 

2 

3from functools import reduce 

4 

5import numpy as np 

6from scipy import stats 

7from statsmodels.base.data import handle_data 

8from statsmodels.tools.data import _is_using_pandas 

9from statsmodels.tools.tools import recipr, nan_dot 

10from statsmodels.stats.contrast import (ContrastResults, WaldTestResults, 

11 t_test_pairwise) 

12from statsmodels.tools.decorators import (cache_readonly, 

13 cached_value, cached_data) 

14import statsmodels.base.wrapper as wrap 

15from statsmodels.tools.numdiff import approx_fprime 

16from statsmodels.tools.sm_exceptions import ValueWarning, \ 

17 HessianInversionWarning 

18from statsmodels.formula import handle_formula_data 

19from statsmodels.base.optimizer import Optimizer 

20 

21 

22_model_params_doc = """Parameters 

23 ---------- 

24 endog : array_like 

25 A 1-d endogenous response variable. The dependent variable. 

26 exog : array_like 

27 A nobs x k array where `nobs` is the number of observations and `k` 

28 is the number of regressors. An intercept is not included by default 

29 and should be added by the user. See 

30 :func:`statsmodels.tools.add_constant`.""" 

31 

32_missing_param_doc = """\ 

33missing : str 

34 Available options are 'none', 'drop', and 'raise'. If 'none', no nan 

35 checking is done. If 'drop', any observations with nans are dropped. 

36 If 'raise', an error is raised. Default is 'none'.""" 

37_extra_param_doc = """ 

38 hasconst : None or bool 

39 Indicates whether the RHS includes a user-supplied constant. If True, 

40 a constant is not checked for and k_constant is set to 1 and all 

41 result statistics are calculated as if a constant is present. If 

42 False, a constant is not checked for and k_constant is set to 0. 

43 **kwargs 

44 Extra arguments that are used to set model properties when using the 

45 formula interface.""" 

46 

47 

48class Model(object): 

49 __doc__ = """ 

50 A (predictive) statistical model. Intended to be subclassed not used. 

51 

52 %(params_doc)s 

53 %(extra_params_doc)s 

54 

55 Attributes 

56 ---------- 

57 exog_names 

58 endog_names 

59 

60 Notes 

61 ----- 

62 `endog` and `exog` are references to any data provided. So if the data is 

63 already stored in numpy arrays and it is changed then `endog` and `exog` 

64 will change as well. 

65 """ % {'params_doc': _model_params_doc, 

66 'extra_params_doc': _missing_param_doc + _extra_param_doc} 

67 

68 # Maximum number of endogenous variables when using a formula 

69 # Default is 1, which is more common. Override in models when needed 

70 # Set to None to skip check 

71 _formula_max_endog = 1 

72 

73 def __init__(self, endog, exog=None, **kwargs): 

74 missing = kwargs.pop('missing', 'none') 

75 hasconst = kwargs.pop('hasconst', None) 

76 self.data = self._handle_data(endog, exog, missing, hasconst, 

77 **kwargs) 

78 self.k_constant = self.data.k_constant 

79 self.exog = self.data.exog 

80 self.endog = self.data.endog 

81 self._data_attr = [] 

82 self._data_attr.extend(['exog', 'endog', 'data.exog', 'data.endog']) 

83 if 'formula' not in kwargs: # will not be able to unpickle without these 

84 self._data_attr.extend(['data.orig_endog', 'data.orig_exog']) 

85 # store keys for extras if we need to recreate model instance 

86 # we do not need 'missing', maybe we need 'hasconst' 

87 self._init_keys = list(kwargs.keys()) 

88 if hasconst is not None: 

89 self._init_keys.append('hasconst') 

90 

91 def _get_init_kwds(self): 

92 """return dictionary with extra keys used in model.__init__ 

93 """ 

94 kwds = dict(((key, getattr(self, key, None)) 

95 for key in self._init_keys)) 

96 

97 return kwds 

98 

99 def _handle_data(self, endog, exog, missing, hasconst, **kwargs): 

100 data = handle_data(endog, exog, missing, hasconst, **kwargs) 

101 # kwargs arrays could have changed, easier to just attach here 

102 for key in kwargs: 

103 if key in ['design_info', 'formula']: # leave attached to data 

104 continue 

105 # pop so we do not start keeping all these twice or references 

106 try: 

107 setattr(self, key, data.__dict__.pop(key)) 

108 except KeyError: # panel already pops keys in data handling 

109 pass 

110 return data 

111 

112 @classmethod 

113 def from_formula(cls, formula, data, subset=None, drop_cols=None, 

114 *args, **kwargs): 

115 """ 

116 Create a Model from a formula and dataframe. 

117 

118 Parameters 

119 ---------- 

120 formula : str or generic Formula object 

121 The formula specifying the model. 

122 data : array_like 

123 The data for the model. See Notes. 

124 subset : array_like 

125 An array-like object of booleans, integers, or index values that 

126 indicate the subset of df to use in the model. Assumes df is a 

127 `pandas.DataFrame`. 

128 drop_cols : array_like 

129 Columns to drop from the design matrix. Cannot be used to 

130 drop terms involving categoricals. 

131 *args 

132 Additional positional argument that are passed to the model. 

133 **kwargs 

134 These are passed to the model with one exception. The 

135 ``eval_env`` keyword is passed to patsy. It can be either a 

136 :class:`patsy:patsy.EvalEnvironment` object or an integer 

137 indicating the depth of the namespace to use. For example, the 

138 default ``eval_env=0`` uses the calling namespace. If you wish 

139 to use a "clean" environment set ``eval_env=-1``. 

140 

141 Returns 

142 ------- 

143 model 

144 The model instance. 

145 

146 Notes 

147 ----- 

148 data must define __getitem__ with the keys in the formula terms 

149 args and kwargs are passed on to the model instantiation. E.g., 

150 a numpy structured or rec array, a dictionary, or a pandas DataFrame. 

151 """ 

152 # TODO: provide a docs template for args/kwargs from child models 

153 # TODO: subset could use syntax. issue #469. 

154 if subset is not None: 

155 data = data.loc[subset] 

156 eval_env = kwargs.pop('eval_env', None) 

157 if eval_env is None: 

158 eval_env = 2 

159 elif eval_env == -1: 

160 from patsy import EvalEnvironment 

161 eval_env = EvalEnvironment({}) 

162 elif isinstance(eval_env, int): 

163 eval_env += 1 # we're going down the stack again 

164 missing = kwargs.get('missing', 'drop') 

165 if missing == 'none': # with patsy it's drop or raise. let's raise. 

166 missing = 'raise' 

167 

168 tmp = handle_formula_data(data, None, formula, depth=eval_env, 

169 missing=missing) 

170 ((endog, exog), missing_idx, design_info) = tmp 

171 max_endog = cls._formula_max_endog 

172 if (max_endog is not None and 

173 endog.ndim > 1 and endog.shape[1] > max_endog): 

174 raise ValueError('endog has evaluated to an array with multiple ' 

175 'columns that has shape {0}. This occurs when ' 

176 'the variable converted to endog is non-numeric' 

177 ' (e.g., bool or str).'.format(endog.shape)) 

178 if drop_cols is not None and len(drop_cols) > 0: 

179 cols = [x for x in exog.columns if x not in drop_cols] 

180 if len(cols) < len(exog.columns): 

181 exog = exog[cols] 

182 cols = list(design_info.term_names) 

183 for col in drop_cols: 

184 try: 

185 cols.remove(col) 

186 except ValueError: 

187 pass # OK if not present 

188 design_info = design_info.subset(cols) 

189 

190 kwargs.update({'missing_idx': missing_idx, 

191 'missing': missing, 

192 'formula': formula, # attach formula for unpckling 

193 'design_info': design_info}) 

194 mod = cls(endog, exog, *args, **kwargs) 

195 mod.formula = formula 

196 

197 # since we got a dataframe, attach the original 

198 mod.data.frame = data 

199 return mod 

200 

201 @property 

202 def endog_names(self): 

203 """ 

204 Names of endogenous variables. 

205 """ 

206 return self.data.ynames 

207 

208 @property 

209 def exog_names(self): 

210 """ 

211 Names of exogenous variables. 

212 """ 

213 return self.data.xnames 

214 

215 def fit(self): 

216 """ 

217 Fit a model to data. 

218 """ 

219 raise NotImplementedError 

220 

221 def predict(self, params, exog=None, *args, **kwargs): 

222 """ 

223 After a model has been fit predict returns the fitted values. 

224 

225 This is a placeholder intended to be overwritten by individual models. 

226 """ 

227 raise NotImplementedError 

228 

229 

230class LikelihoodModel(Model): 

231 """ 

232 Likelihood model is a subclass of Model. 

233 """ 

234 

235 def __init__(self, endog, exog=None, **kwargs): 

236 super(LikelihoodModel, self).__init__(endog, exog, **kwargs) 

237 self.initialize() 

238 

239 def initialize(self): 

240 """ 

241 Initialize (possibly re-initialize) a Model instance. 

242 

243 For example, if the the design matrix of a linear model changes then 

244 initialized can be used to recompute values using the modified design 

245 matrix. 

246 """ 

247 pass 

248 

249 # TODO: if the intent is to re-initialize the model with new data then this 

250 # method needs to take inputs... 

251 

252 def loglike(self, params): 

253 """ 

254 Log-likelihood of model. 

255 

256 Parameters 

257 ---------- 

258 params : ndarray 

259 The model parameters used to compute the log-likelihood. 

260 

261 Notes 

262 ----- 

263 Must be overridden by subclasses. 

264 """ 

265 raise NotImplementedError 

266 

267 def score(self, params): 

268 """ 

269 Score vector of model. 

270 

271 The gradient of logL with respect to each parameter. 

272 

273 Parameters 

274 ---------- 

275 params : ndarray 

276 The parameters to use when evaluating the Hessian. 

277 

278 Returns 

279 ------- 

280 ndarray 

281 The score vector evaluated at the parameters. 

282 """ 

283 raise NotImplementedError 

284 

285 def information(self, params): 

286 """ 

287 Fisher information matrix of model. 

288 

289 Returns -1 * Hessian of the log-likelihood evaluated at params. 

290 

291 Parameters 

292 ---------- 

293 params : ndarray 

294 The model parameters. 

295 """ 

296 raise NotImplementedError 

297 

298 def hessian(self, params): 

299 """ 

300 The Hessian matrix of the model. 

301 

302 Parameters 

303 ---------- 

304 params : ndarray 

305 The parameters to use when evaluating the Hessian. 

306 

307 Returns 

308 ------- 

309 ndarray 

310 The hessian evaluated at the parameters. 

311 """ 

312 raise NotImplementedError 

313 

314 def fit(self, start_params=None, method='newton', maxiter=100, 

315 full_output=True, disp=True, fargs=(), callback=None, retall=False, 

316 skip_hessian=False, **kwargs): 

317 """ 

318 Fit method for likelihood based models 

319 

320 Parameters 

321 ---------- 

322 start_params : array_like, optional 

323 Initial guess of the solution for the loglikelihood maximization. 

324 The default is an array of zeros. 

325 method : str, optional 

326 The `method` determines which solver from `scipy.optimize` 

327 is used, and it can be chosen from among the following strings: 

328 

329 - 'newton' for Newton-Raphson, 'nm' for Nelder-Mead 

330 - 'bfgs' for Broyden-Fletcher-Goldfarb-Shanno (BFGS) 

331 - 'lbfgs' for limited-memory BFGS with optional box constraints 

332 - 'powell' for modified Powell's method 

333 - 'cg' for conjugate gradient 

334 - 'ncg' for Newton-conjugate gradient 

335 - 'basinhopping' for global basin-hopping solver 

336 - 'minimize' for generic wrapper of scipy minimize (BFGS by default) 

337 

338 The explicit arguments in `fit` are passed to the solver, 

339 with the exception of the basin-hopping solver. Each 

340 solver has several optional arguments that are not the same across 

341 solvers. See the notes section below (or scipy.optimize) for the 

342 available arguments and for the list of explicit arguments that the 

343 basin-hopping solver supports. 

344 maxiter : int, optional 

345 The maximum number of iterations to perform. 

346 full_output : bool, optional 

347 Set to True to have all available output in the Results object's 

348 mle_retvals attribute. The output is dependent on the solver. 

349 See LikelihoodModelResults notes section for more information. 

350 disp : bool, optional 

351 Set to True to print convergence messages. 

352 fargs : tuple, optional 

353 Extra arguments passed to the likelihood function, i.e., 

354 loglike(x,*args) 

355 callback : callable callback(xk), optional 

356 Called after each iteration, as callback(xk), where xk is the 

357 current parameter vector. 

358 retall : bool, optional 

359 Set to True to return list of solutions at each iteration. 

360 Available in Results object's mle_retvals attribute. 

361 skip_hessian : bool, optional 

362 If False (default), then the negative inverse hessian is calculated 

363 after the optimization. If True, then the hessian will not be 

364 calculated. However, it will be available in methods that use the 

365 hessian in the optimization (currently only with `"newton"`). 

366 kwargs : keywords 

367 All kwargs are passed to the chosen solver with one exception. The 

368 following keyword controls what happens after the fit:: 

369 

370 warn_convergence : bool, optional 

371 If True, checks the model for the converged flag. If the 

372 converged flag is False, a ConvergenceWarning is issued. 

373 

374 Notes 

375 ----- 

376 The 'basinhopping' solver ignores `maxiter`, `retall`, `full_output` 

377 explicit arguments. 

378 

379 Optional arguments for solvers (see returned Results.mle_settings):: 

380 

381 'newton' 

382 tol : float 

383 Relative error in params acceptable for convergence. 

384 'nm' -- Nelder Mead 

385 xtol : float 

386 Relative error in params acceptable for convergence 

387 ftol : float 

388 Relative error in loglike(params) acceptable for 

389 convergence 

390 maxfun : int 

391 Maximum number of function evaluations to make. 

392 'bfgs' 

393 gtol : float 

394 Stop when norm of gradient is less than gtol. 

395 norm : float 

396 Order of norm (np.Inf is max, -np.Inf is min) 

397 epsilon 

398 If fprime is approximated, use this value for the step 

399 size. Only relevant if LikelihoodModel.score is None. 

400 'lbfgs' 

401 m : int 

402 This many terms are used for the Hessian approximation. 

403 factr : float 

404 A stop condition that is a variant of relative error. 

405 pgtol : float 

406 A stop condition that uses the projected gradient. 

407 epsilon 

408 If fprime is approximated, use this value for the step 

409 size. Only relevant if LikelihoodModel.score is None. 

410 maxfun : int 

411 Maximum number of function evaluations to make. 

412 bounds : sequence 

413 (min, max) pairs for each element in x, 

414 defining the bounds on that parameter. 

415 Use None for one of min or max when there is no bound 

416 in that direction. 

417 'cg' 

418 gtol : float 

419 Stop when norm of gradient is less than gtol. 

420 norm : float 

421 Order of norm (np.Inf is max, -np.Inf is min) 

422 epsilon : float 

423 If fprime is approximated, use this value for the step 

424 size. Can be scalar or vector. Only relevant if 

425 Likelihoodmodel.score is None. 

426 'ncg' 

427 fhess_p : callable f'(x,*args) 

428 Function which computes the Hessian of f times an arbitrary 

429 vector, p. Should only be supplied if 

430 LikelihoodModel.hessian is None. 

431 avextol : float 

432 Stop when the average relative error in the minimizer 

433 falls below this amount. 

434 epsilon : float or ndarray 

435 If fhess is approximated, use this value for the step size. 

436 Only relevant if Likelihoodmodel.hessian is None. 

437 'powell' 

438 xtol : float 

439 Line-search error tolerance 

440 ftol : float 

441 Relative error in loglike(params) for acceptable for 

442 convergence. 

443 maxfun : int 

444 Maximum number of function evaluations to make. 

445 start_direc : ndarray 

446 Initial direction set. 

447 'basinhopping' 

448 niter : int 

449 The number of basin hopping iterations. 

450 niter_success : int 

451 Stop the run if the global minimum candidate remains the 

452 same for this number of iterations. 

453 T : float 

454 The "temperature" parameter for the accept or reject 

455 criterion. Higher "temperatures" mean that larger jumps 

456 in function value will be accepted. For best results 

457 `T` should be comparable to the separation (in function 

458 value) between local minima. 

459 stepsize : float 

460 Initial step size for use in the random displacement. 

461 interval : int 

462 The interval for how often to update the `stepsize`. 

463 minimizer : dict 

464 Extra keyword arguments to be passed to the minimizer 

465 `scipy.optimize.minimize()`, for example 'method' - the 

466 minimization method (e.g. 'L-BFGS-B'), or 'tol' - the 

467 tolerance for termination. Other arguments are mapped from 

468 explicit argument of `fit`: 

469 - `args` <- `fargs` 

470 - `jac` <- `score` 

471 - `hess` <- `hess` 

472 'minimize' 

473 min_method : str, optional 

474 Name of minimization method to use. 

475 Any method specific arguments can be passed directly. 

476 For a list of methods and their arguments, see 

477 documentation of `scipy.optimize.minimize`. 

478 If no method is specified, then BFGS is used. 

479 """ 

480 Hinv = None # JP error if full_output=0, Hinv not defined 

481 

482 if start_params is None: 

483 if hasattr(self, 'start_params'): 

484 start_params = self.start_params 

485 elif self.exog is not None: 

486 # fails for shape (K,)? 

487 start_params = [0] * self.exog.shape[1] 

488 else: 

489 raise ValueError("If exog is None, then start_params should " 

490 "be specified") 

491 

492 # TODO: separate args from nonarg taking score and hessian, ie., 

493 # user-supplied and numerically evaluated estimate frprime does not take 

494 # args in most (any?) of the optimize function 

495 

496 nobs = self.endog.shape[0] 

497 # f = lambda params, *args: -self.loglike(params, *args) / nobs 

498 

499 def f(params, *args): 

500 return -self.loglike(params, *args) / nobs 

501 

502 if method == 'newton': 

503 # TODO: why are score and hess positive? 

504 def score(params, *args): 

505 return self.score(params, *args) / nobs 

506 

507 def hess(params, *args): 

508 return self.hessian(params, *args) / nobs 

509 else: 

510 def score(params, *args): 

511 return -self.score(params, *args) / nobs 

512 

513 def hess(params, *args): 

514 return -self.hessian(params, *args) / nobs 

515 

516 warn_convergence = kwargs.pop('warn_convergence', True) 

517 optimizer = Optimizer() 

518 xopt, retvals, optim_settings = optimizer._fit(f, score, start_params, 

519 fargs, kwargs, 

520 hessian=hess, 

521 method=method, 

522 disp=disp, 

523 maxiter=maxiter, 

524 callback=callback, 

525 retall=retall, 

526 full_output=full_output) 

527 

528 # NOTE: this is for fit_regularized and should be generalized 

529 cov_params_func = kwargs.setdefault('cov_params_func', None) 

530 if cov_params_func: 

531 Hinv = cov_params_func(self, xopt, retvals) 

532 elif method == 'newton' and full_output: 

533 Hinv = np.linalg.inv(-retvals['Hessian']) / nobs 

534 elif not skip_hessian: 

535 H = -1 * self.hessian(xopt) 

536 invertible = False 

537 if np.all(np.isfinite(H)): 

538 eigvals, eigvecs = np.linalg.eigh(H) 

539 if np.min(eigvals) > 0: 

540 invertible = True 

541 

542 if invertible: 

543 Hinv = eigvecs.dot(np.diag(1.0 / eigvals)).dot(eigvecs.T) 

544 Hinv = np.asfortranarray((Hinv + Hinv.T) / 2.0) 

545 else: 

546 from warnings import warn 

547 warn('Inverting hessian failed, no bse or cov_params ' 

548 'available', HessianInversionWarning) 

549 Hinv = None 

550 

551 if 'cov_type' in kwargs: 

552 cov_kwds = kwargs.get('cov_kwds', {}) 

553 kwds = {'cov_type': kwargs['cov_type'], 'cov_kwds': cov_kwds} 

554 else: 

555 kwds = {} 

556 if 'use_t' in kwargs: 

557 kwds['use_t'] = kwargs['use_t'] 

558 # TODO: add Hessian approximation and change the above if needed 

559 mlefit = LikelihoodModelResults(self, xopt, Hinv, scale=1., **kwds) 

560 

561 # TODO: hardcode scale? 

562 mlefit.mle_retvals = retvals 

563 if isinstance(retvals, dict): 

564 if warn_convergence and not retvals['converged']: 

565 from warnings import warn 

566 from statsmodels.tools.sm_exceptions import ConvergenceWarning 

567 warn("Maximum Likelihood optimization failed to converge. " 

568 "Check mle_retvals", ConvergenceWarning) 

569 

570 mlefit.mle_settings = optim_settings 

571 return mlefit 

572 

573 def _fit_zeros(self, keep_index=None, start_params=None, 

574 return_auxiliary=False, k_params=None, **fit_kwds): 

575 """experimental, fit the model subject to zero constraints 

576 

577 Intended for internal use cases until we know what we need. 

578 API will need to change to handle models with two exog. 

579 This is not yet supported by all model subclasses. 

580 

581 This is essentially a simplified version of `fit_constrained`, and 

582 does not need to use `offset`. 

583 

584 The estimation creates a new model with transformed design matrix, 

585 exog, and converts the results back to the original parameterization. 

586 

587 Some subclasses could use a more efficient calculation than using a 

588 new model. 

589 

590 Parameters 

591 ---------- 

592 keep_index : array_like (int or bool) or slice 

593 variables that should be dropped. 

594 start_params : None or array_like 

595 starting values for the optimization. `start_params` needs to be 

596 given in the original parameter space and are internally 

597 transformed. 

598 k_params : int or None 

599 If None, then we try to infer from start_params or model. 

600 **fit_kwds : keyword arguments 

601 fit_kwds are used in the optimization of the transformed model. 

602 

603 Returns 

604 ------- 

605 results : Results instance 

606 """ 

607 # we need to append index of extra params to keep_index as in 

608 # NegativeBinomial 

609 if hasattr(self, 'k_extra') and self.k_extra > 0: 

610 # we cannot change the original, TODO: should we add keep_index_params? 

611 keep_index = np.array(keep_index, copy=True) 

612 k = self.exog.shape[1] 

613 extra_index = np.arange(k, k + self.k_extra) 

614 keep_index_p = np.concatenate((keep_index, extra_index)) 

615 else: 

616 keep_index_p = keep_index 

617 

618 # not all models support start_params, drop if None, hide them in fit_kwds 

619 if start_params is not None: 

620 fit_kwds['start_params'] = start_params[keep_index_p] 

621 k_params = len(start_params) 

622 # ignore k_params in this case, or verify consisteny? 

623 

624 # build auxiliary model and fit 

625 init_kwds = self._get_init_kwds() 

626 mod_constr = self.__class__(self.endog, self.exog[:, keep_index], 

627 **init_kwds) 

628 res_constr = mod_constr.fit(**fit_kwds) 

629 # switch name, only need keep_index for params below 

630 keep_index = keep_index_p 

631 

632 if k_params is None: 

633 k_params = self.exog.shape[1] 

634 k_params += getattr(self, 'k_extra', 0) 

635 

636 params_full = np.zeros(k_params) 

637 params_full[keep_index] = res_constr.params 

638 

639 # create dummy results Instance, TODO: wire up properly 

640 # TODO: this could be moved into separate private method if needed 

641 # discrete L1 fit_regularized doens't reestimate AFAICS 

642 # RLM does not have method, disp nor warn_convergence keywords 

643 # OLS, WLS swallows extra kwds with **kwargs, but does not have method='nm' 

644 try: 

645 # Note: addding full_output=False causes exceptions 

646 res = self.fit(maxiter=0, disp=0, method='nm', skip_hessian=True, 

647 warn_convergence=False, start_params=params_full) 

648 # we get a wrapper back 

649 except (TypeError, ValueError): 

650 res = self.fit() 

651 

652 # Warning: make sure we are not just changing the wrapper instead of 

653 # results #2400 

654 # TODO: do we need to change res._results.scale in some models? 

655 if hasattr(res_constr.model, 'scale'): 

656 # Note: res.model is self 

657 # GLM problem, see #2399, 

658 # TODO: remove from model if not needed anymore 

659 res.model.scale = res._results.scale = res_constr.model.scale 

660 

661 if hasattr(res_constr, 'mle_retvals'): 

662 res._results.mle_retvals = res_constr.mle_retvals 

663 # not available for not scipy optimization, e.g. glm irls 

664 # TODO: what retvals should be required? 

665 # res.mle_retvals['fcall'] = res_constr.mle_retvals.get('fcall', np.nan) 

666 # res.mle_retvals['iterations'] = res_constr.mle_retvals.get( 

667 # 'iterations', np.nan) 

668 # res.mle_retvals['converged'] = res_constr.mle_retvals['converged'] 

669 # overwrite all mle_settings 

670 if hasattr(res_constr, 'mle_settings'): 

671 res._results.mle_settings = res_constr.mle_settings 

672 

673 res._results.params = params_full 

674 if (not hasattr(res._results, 'normalized_cov_params') or 

675 res._results.normalized_cov_params is None): 

676 res._results.normalized_cov_params = np.zeros((k_params, k_params)) 

677 else: 

678 res._results.normalized_cov_params[...] = 0 

679 

680 # fancy indexing requires integer array 

681 keep_index = np.array(keep_index) 

682 res._results.normalized_cov_params[keep_index[:, None], keep_index] = \ 

683 res_constr.normalized_cov_params 

684 k_constr = res_constr.df_resid - res._results.df_resid 

685 if hasattr(res_constr, 'cov_params_default'): 

686 res._results.cov_params_default = np.zeros((k_params, k_params)) 

687 res._results.cov_params_default[keep_index[:, None], keep_index] = \ 

688 res_constr.cov_params_default 

689 if hasattr(res_constr, 'cov_type'): 

690 res._results.cov_type = res_constr.cov_type 

691 res._results.cov_kwds = res_constr.cov_kwds 

692 

693 res._results.keep_index = keep_index 

694 res._results.df_resid = res_constr.df_resid 

695 res._results.df_model = res_constr.df_model 

696 

697 res._results.k_constr = k_constr 

698 res._results.results_constrained = res_constr 

699 

700 # special temporary workaround for RLM 

701 # need to be able to override robust covariances 

702 if hasattr(res.model, 'M'): 

703 del res._results._cache['resid'] 

704 del res._results._cache['fittedvalues'] 

705 del res._results._cache['sresid'] 

706 cov = res._results._cache['bcov_scaled'] 

707 # inplace adjustment 

708 cov[...] = 0 

709 cov[keep_index[:, None], keep_index] = res_constr.bcov_scaled 

710 res._results.cov_params_default = cov 

711 

712 return res 

713 

714 def _fit_collinear(self, atol=1e-14, rtol=1e-13, **kwds): 

715 """experimental, fit of the model without collinear variables 

716 

717 This currently uses QR to drop variables based on the given 

718 sequence. 

719 Options will be added in future, when the supporting functions 

720 to identify collinear variables become available. 

721 """ 

722 

723 # ------ copied from PR #2380 remove when merged 

724 x = self.exog 

725 tol = atol + rtol * x.var(0) 

726 r = np.linalg.qr(x, mode='r') 

727 mask = np.abs(r.diagonal()) < np.sqrt(tol) 

728 # TODO add to results instance 

729 # idx_collinear = np.where(mask)[0] 

730 idx_keep = np.where(~mask)[0] 

731 return self._fit_zeros(keep_index=idx_keep, **kwds) 

732 

733 

734# TODO: the below is unfinished 

735class GenericLikelihoodModel(LikelihoodModel): 

736 """ 

737 Allows the fitting of any likelihood function via maximum likelihood. 

738 

739 A subclass needs to specify at least the log-likelihood 

740 If the log-likelihood is specified for each observation, then results that 

741 require the Jacobian will be available. (The other case is not tested yet.) 

742 

743 Notes 

744 ----- 

745 Optimization methods that require only a likelihood function are 'nm' and 

746 'powell' 

747 

748 Optimization methods that require a likelihood function and a 

749 score/gradient are 'bfgs', 'cg', and 'ncg'. A function to compute the 

750 Hessian is optional for 'ncg'. 

751 

752 Optimization method that require a likelihood function, a score/gradient, 

753 and a Hessian is 'newton' 

754 

755 If they are not overwritten by a subclass, then numerical gradient, 

756 Jacobian and Hessian of the log-likelihood are calculated by numerical 

757 forward differentiation. This might results in some cases in precision 

758 problems, and the Hessian might not be positive definite. Even if the 

759 Hessian is not positive definite the covariance matrix of the parameter 

760 estimates based on the outer product of the Jacobian might still be valid. 

761 

762 

763 Examples 

764 -------- 

765 see also subclasses in directory miscmodels 

766 

767 import statsmodels.api as sm 

768 data = sm.datasets.spector.load(as_pandas=False) 

769 data.exog = sm.add_constant(data.exog) 

770 # in this dir 

771 from model import GenericLikelihoodModel 

772 probit_mod = sm.Probit(data.endog, data.exog) 

773 probit_res = probit_mod.fit() 

774 loglike = probit_mod.loglike 

775 score = probit_mod.score 

776 mod = GenericLikelihoodModel(data.endog, data.exog, loglike, score) 

777 res = mod.fit(method="nm", maxiter = 500) 

778 import numpy as np 

779 np.allclose(res.params, probit_res.params) 

780 """ 

781 def __init__(self, endog, exog=None, loglike=None, score=None, 

782 hessian=None, missing='none', extra_params_names=None, 

783 **kwds): 

784 # let them be none in case user wants to use inheritance 

785 if loglike is not None: 

786 self.loglike = loglike 

787 if score is not None: 

788 self.score = score 

789 if hessian is not None: 

790 self.hessian = hessian 

791 

792 self.__dict__.update(kwds) 

793 

794 # TODO: data structures? 

795 

796 # TODO temporary solution, force approx normal 

797 # self.df_model = 9999 

798 # somewhere: CacheWriteWarning: 'df_model' cannot be overwritten 

799 super(GenericLikelihoodModel, self).__init__(endog, exog, 

800 missing=missing) 

801 

802 # this will not work for ru2nmnl, maybe np.ndim of a dict? 

803 if exog is not None: 

804 self.nparams = (exog.shape[1] if np.ndim(exog) == 2 else 1) 

805 

806 if extra_params_names is not None: 

807 self._set_extra_params_names(extra_params_names) 

808 

809 def _set_extra_params_names(self, extra_params_names): 

810 # check param_names 

811 if extra_params_names is not None: 

812 if self.exog is not None: 

813 self.exog_names.extend(extra_params_names) 

814 else: 

815 self.data.xnames = extra_params_names 

816 

817 self.nparams = len(self.exog_names) 

818 

819 # this is redundant and not used when subclassing 

820 def initialize(self): 

821 """ 

822 Initialize (possibly re-initialize) a Model instance. For 

823 instance, the design matrix of a linear model may change 

824 and some things must be recomputed. 

825 """ 

826 if not self.score: # right now score is not optional 

827 self.score = lambda x: approx_fprime(x, self.loglike) 

828 if not self.hessian: 

829 pass 

830 else: # can use approx_hess_p if we have a gradient 

831 if not self.hessian: 

832 pass 

833 # Initialize is called by 

834 # statsmodels.model.LikelihoodModel.__init__ 

835 # and should contain any preprocessing that needs to be done for a model 

836 if self.exog is not None: 

837 # assume constant 

838 er = np.linalg.matrix_rank(self.exog) 

839 self.df_model = float(er - 1) 

840 self.df_resid = float(self.exog.shape[0] - er) 

841 else: 

842 self.df_model = np.nan 

843 self.df_resid = np.nan 

844 super(GenericLikelihoodModel, self).initialize() 

845 

846 def expandparams(self, params): 

847 """ 

848 expand to full parameter array when some parameters are fixed 

849 

850 Parameters 

851 ---------- 

852 params : ndarray 

853 reduced parameter array 

854 

855 Returns 

856 ------- 

857 paramsfull : ndarray 

858 expanded parameter array where fixed parameters are included 

859 

860 Notes 

861 ----- 

862 Calling this requires that self.fixed_params and self.fixed_paramsmask 

863 are defined. 

864 

865 *developer notes:* 

866 

867 This can be used in the log-likelihood to ... 

868 

869 this could also be replaced by a more general parameter 

870 transformation. 

871 """ 

872 paramsfull = self.fixed_params.copy() 

873 paramsfull[self.fixed_paramsmask] = params 

874 return paramsfull 

875 

876 def reduceparams(self, params): 

877 """Reduce parameters""" 

878 return params[self.fixed_paramsmask] 

879 

880 def loglike(self, params): 

881 """Log-likelihood of model at params""" 

882 return self.loglikeobs(params).sum(0) 

883 

884 def nloglike(self, params): 

885 """Negative log-likelihood of model at params""" 

886 return -self.loglikeobs(params).sum(0) 

887 

888 def loglikeobs(self, params): 

889 """Log-likelihood of individual observations at params""" 

890 return -self.nloglikeobs(params) 

891 

892 def score(self, params): 

893 """ 

894 Gradient of log-likelihood evaluated at params 

895 """ 

896 kwds = {} 

897 kwds.setdefault('centered', True) 

898 return approx_fprime(params, self.loglike, **kwds).ravel() 

899 

900 def score_obs(self, params, **kwds): 

901 """ 

902 Jacobian/Gradient of log-likelihood evaluated at params for each 

903 observation. 

904 """ 

905 # kwds.setdefault('epsilon', 1e-4) 

906 kwds.setdefault('centered', True) 

907 return approx_fprime(params, self.loglikeobs, **kwds) 

908 

909 def hessian(self, params): 

910 """ 

911 Hessian of log-likelihood evaluated at params 

912 """ 

913 from statsmodels.tools.numdiff import approx_hess 

914 # need options for hess (epsilon) 

915 return approx_hess(params, self.loglike) 

916 

917 def hessian_factor(self, params, scale=None, observed=True): 

918 """Weights for calculating Hessian 

919 

920 Parameters 

921 ---------- 

922 params : ndarray 

923 parameter at which Hessian is evaluated 

924 scale : None or float 

925 If scale is None, then the default scale will be calculated. 

926 Default scale is defined by `self.scaletype` and set in fit. 

927 If scale is not None, then it is used as a fixed scale. 

928 observed : bool 

929 If True, then the observed Hessian is returned. If false then the 

930 expected information matrix is returned. 

931 

932 Returns 

933 ------- 

934 hessian_factor : ndarray, 1d 

935 A 1d weight vector used in the calculation of the Hessian. 

936 The hessian is obtained by `(exog.T * hessian_factor).dot(exog)` 

937 """ 

938 

939 raise NotImplementedError 

940 

941 def fit(self, start_params=None, method='nm', maxiter=500, full_output=1, 

942 disp=1, callback=None, retall=0, **kwargs): 

943 """ 

944 Fit the model using maximum likelihood. 

945 

946 See LikelihoodModel.fit for more information. 

947 """ 

948 if start_params is None: 

949 if hasattr(self, 'start_params'): 

950 start_params = self.start_params 

951 else: 

952 start_params = 0.1 * np.ones(self.nparams) 

953 

954 fit_method = super(GenericLikelihoodModel, self).fit 

955 mlefit = fit_method(start_params=start_params, 

956 method=method, maxiter=maxiter, 

957 full_output=full_output, 

958 disp=disp, callback=callback, **kwargs) 

959 genericmlefit = GenericLikelihoodModelResults(self, mlefit) 

960 

961 # amend param names 

962 exog_names = [] if (self.exog_names is None) else self.exog_names 

963 k_miss = len(exog_names) - len(mlefit.params) 

964 if not k_miss == 0: 

965 if k_miss < 0: 

966 self._set_extra_params_names(['par%d' % i 

967 for i in range(-k_miss)]) 

968 else: 

969 # I do not want to raise after we have already fit() 

970 import warnings 

971 warnings.warn('more exog_names than parameters', ValueWarning) 

972 

973 return genericmlefit 

974 

975 

976class Results(object): 

977 """ 

978 Class to contain model results 

979 

980 Parameters 

981 ---------- 

982 model : class instance 

983 the previously specified model instance 

984 params : ndarray 

985 parameter estimates from the fit model 

986 """ 

987 def __init__(self, model, params, **kwd): 

988 self.__dict__.update(kwd) 

989 self.initialize(model, params, **kwd) 

990 self._data_attr = [] 

991 

992 def initialize(self, model, params, **kwargs): 

993 """ 

994 Initialize (possibly re-initialize) a Results instance. 

995 

996 Parameters 

997 ---------- 

998 model : Model 

999 The model instance. 

1000 params : ndarray 

1001 The model parameters. 

1002 **kwargs 

1003 Any additional keyword arguments required to initialize the model. 

1004 """ 

1005 self.params = params 

1006 self.model = model 

1007 if hasattr(model, 'k_constant'): 

1008 self.k_constant = model.k_constant 

1009 

1010 def predict(self, exog=None, transform=True, *args, **kwargs): 

1011 """ 

1012 Call self.model.predict with self.params as the first argument. 

1013 

1014 Parameters 

1015 ---------- 

1016 exog : array_like, optional 

1017 The values for which you want to predict. see Notes below. 

1018 transform : bool, optional 

1019 If the model was fit via a formula, do you want to pass 

1020 exog through the formula. Default is True. E.g., if you fit 

1021 a model y ~ log(x1) + log(x2), and transform is True, then 

1022 you can pass a data structure that contains x1 and x2 in 

1023 their original form. Otherwise, you'd need to log the data 

1024 first. 

1025 *args 

1026 Additional arguments to pass to the model, see the 

1027 predict method of the model for the details. 

1028 **kwargs 

1029 Additional keywords arguments to pass to the model, see the 

1030 predict method of the model for the details. 

1031 

1032 Returns 

1033 ------- 

1034 array_like 

1035 See self.model.predict. 

1036 

1037 Notes 

1038 ----- 

1039 The types of exog that are supported depends on whether a formula 

1040 was used in the specification of the model. 

1041 

1042 If a formula was used, then exog is processed in the same way as 

1043 the original data. This transformation needs to have key access to the 

1044 same variable names, and can be a pandas DataFrame or a dict like 

1045 object that contains numpy arrays. 

1046 

1047 If no formula was used, then the provided exog needs to have the 

1048 same number of columns as the original exog in the model. No 

1049 transformation of the data is performed except converting it to 

1050 a numpy array. 

1051 

1052 Row indices as in pandas data frames are supported, and added to the 

1053 returned prediction. 

1054 """ 

1055 import pandas as pd 

1056 

1057 is_pandas = _is_using_pandas(exog, None) 

1058 

1059 exog_index = exog.index if is_pandas else None 

1060 

1061 if transform and hasattr(self.model, 'formula') and (exog is not None): 

1062 design_info = self.model.data.design_info 

1063 from patsy import dmatrix 

1064 if isinstance(exog, pd.Series): 

1065 # we are guessing whether it should be column or row 

1066 if (hasattr(exog, 'name') and isinstance(exog.name, str) and 

1067 exog.name in design_info.describe()): 

1068 # assume we need one column 

1069 exog = pd.DataFrame(exog) 

1070 else: 

1071 # assume we need a row 

1072 exog = pd.DataFrame(exog).T 

1073 orig_exog_len = len(exog) 

1074 is_dict = isinstance(exog, dict) 

1075 try: 

1076 exog = dmatrix(design_info, exog, return_type="dataframe") 

1077 except Exception as exc: 

1078 msg = ('predict requires that you use a DataFrame when ' 

1079 'predicting from a model\nthat was created using the ' 

1080 'formula api.' 

1081 '\n\nThe original error message returned by patsy is:\n' 

1082 '{0}'.format(str(str(exc)))) 

1083 raise exc.__class__(msg) 

1084 if orig_exog_len > len(exog) and not is_dict: 

1085 import warnings 

1086 if exog_index is None: 

1087 warnings.warn('nan values have been dropped', ValueWarning) 

1088 else: 

1089 exog = exog.reindex(exog_index) 

1090 exog_index = exog.index 

1091 

1092 if exog is not None: 

1093 exog = np.asarray(exog) 

1094 if exog.ndim == 1 and (self.model.exog.ndim == 1 or 

1095 self.model.exog.shape[1] == 1): 

1096 exog = exog[:, None] 

1097 exog = np.atleast_2d(exog) # needed in count model shape[1] 

1098 

1099 predict_results = self.model.predict(self.params, exog, *args, 

1100 **kwargs) 

1101 

1102 if exog_index is not None and not hasattr(predict_results, 

1103 'predicted_values'): 

1104 if predict_results.ndim == 1: 

1105 return pd.Series(predict_results, index=exog_index) 

1106 else: 

1107 return pd.DataFrame(predict_results, index=exog_index) 

1108 else: 

1109 return predict_results 

1110 

1111 def summary(self): 

1112 """ 

1113 Summary 

1114 

1115 Not implemented 

1116 """ 

1117 raise NotImplementedError 

1118 

1119 

1120# TODO: public method? 

1121class LikelihoodModelResults(Results): 

1122 """ 

1123 Class to contain results from likelihood models 

1124 

1125 Parameters 

1126 ---------- 

1127 model : LikelihoodModel instance or subclass instance 

1128 LikelihoodModelResults holds a reference to the model that is fit. 

1129 params : 1d array_like 

1130 parameter estimates from estimated model 

1131 normalized_cov_params : 2d array 

1132 Normalized (before scaling) covariance of params. (dot(X.T,X))**-1 

1133 scale : float 

1134 For (some subset of models) scale will typically be the 

1135 mean square error from the estimated model (sigma^2) 

1136 

1137 Attributes 

1138 ---------- 

1139 mle_retvals : dict 

1140 Contains the values returned from the chosen optimization method if 

1141 full_output is True during the fit. Available only if the model 

1142 is fit by maximum likelihood. See notes below for the output from 

1143 the different methods. 

1144 mle_settings : dict 

1145 Contains the arguments passed to the chosen optimization method. 

1146 Available if the model is fit by maximum likelihood. See 

1147 LikelihoodModel.fit for more information. 

1148 model : model instance 

1149 LikelihoodResults contains a reference to the model that is fit. 

1150 params : ndarray 

1151 The parameters estimated for the model. 

1152 scale : float 

1153 The scaling factor of the model given during instantiation. 

1154 tvalues : ndarray 

1155 The t-values of the standard errors. 

1156 

1157 

1158 Notes 

1159 ----- 

1160 The covariance of params is given by scale times normalized_cov_params. 

1161 

1162 Return values by solver if full_output is True during fit: 

1163 

1164 'newton' 

1165 fopt : float 

1166 The value of the (negative) loglikelihood at its 

1167 minimum. 

1168 iterations : int 

1169 Number of iterations performed. 

1170 score : ndarray 

1171 The score vector at the optimum. 

1172 Hessian : ndarray 

1173 The Hessian at the optimum. 

1174 warnflag : int 

1175 1 if maxiter is exceeded. 0 if successful convergence. 

1176 converged : bool 

1177 True: converged. False: did not converge. 

1178 allvecs : list 

1179 List of solutions at each iteration. 

1180 'nm' 

1181 fopt : float 

1182 The value of the (negative) loglikelihood at its 

1183 minimum. 

1184 iterations : int 

1185 Number of iterations performed. 

1186 warnflag : int 

1187 1: Maximum number of function evaluations made. 

1188 2: Maximum number of iterations reached. 

1189 converged : bool 

1190 True: converged. False: did not converge. 

1191 allvecs : list 

1192 List of solutions at each iteration. 

1193 'bfgs' 

1194 fopt : float 

1195 Value of the (negative) loglikelihood at its minimum. 

1196 gopt : float 

1197 Value of gradient at minimum, which should be near 0. 

1198 Hinv : ndarray 

1199 value of the inverse Hessian matrix at minimum. Note 

1200 that this is just an approximation and will often be 

1201 different from the value of the analytic Hessian. 

1202 fcalls : int 

1203 Number of calls to loglike. 

1204 gcalls : int 

1205 Number of calls to gradient/score. 

1206 warnflag : int 

1207 1: Maximum number of iterations exceeded. 2: Gradient 

1208 and/or function calls are not changing. 

1209 converged : bool 

1210 True: converged. False: did not converge. 

1211 allvecs : list 

1212 Results at each iteration. 

1213 'lbfgs' 

1214 fopt : float 

1215 Value of the (negative) loglikelihood at its minimum. 

1216 gopt : float 

1217 Value of gradient at minimum, which should be near 0. 

1218 fcalls : int 

1219 Number of calls to loglike. 

1220 warnflag : int 

1221 Warning flag: 

1222 

1223 - 0 if converged 

1224 - 1 if too many function evaluations or too many iterations 

1225 - 2 if stopped for another reason 

1226 

1227 converged : bool 

1228 True: converged. False: did not converge. 

1229 'powell' 

1230 fopt : float 

1231 Value of the (negative) loglikelihood at its minimum. 

1232 direc : ndarray 

1233 Current direction set. 

1234 iterations : int 

1235 Number of iterations performed. 

1236 fcalls : int 

1237 Number of calls to loglike. 

1238 warnflag : int 

1239 1: Maximum number of function evaluations. 2: Maximum number 

1240 of iterations. 

1241 converged : bool 

1242 True : converged. False: did not converge. 

1243 allvecs : list 

1244 Results at each iteration. 

1245 'cg' 

1246 fopt : float 

1247 Value of the (negative) loglikelihood at its minimum. 

1248 fcalls : int 

1249 Number of calls to loglike. 

1250 gcalls : int 

1251 Number of calls to gradient/score. 

1252 warnflag : int 

1253 1: Maximum number of iterations exceeded. 2: Gradient and/ 

1254 or function calls not changing. 

1255 converged : bool 

1256 True: converged. False: did not converge. 

1257 allvecs : list 

1258 Results at each iteration. 

1259 'ncg' 

1260 fopt : float 

1261 Value of the (negative) loglikelihood at its minimum. 

1262 fcalls : int 

1263 Number of calls to loglike. 

1264 gcalls : int 

1265 Number of calls to gradient/score. 

1266 hcalls : int 

1267 Number of calls to hessian. 

1268 warnflag : int 

1269 1: Maximum number of iterations exceeded. 

1270 converged : bool 

1271 True: converged. False: did not converge. 

1272 allvecs : list 

1273 Results at each iteration. 

1274 """ 

1275 

1276 # by default we use normal distribution 

1277 # can be overwritten by instances or subclasses 

1278 

1279 def __init__(self, model, params, normalized_cov_params=None, scale=1., 

1280 **kwargs): 

1281 super(LikelihoodModelResults, self).__init__(model, params) 

1282 self.normalized_cov_params = normalized_cov_params 

1283 self.scale = scale 

1284 self._use_t = False 

1285 # robust covariance 

1286 # We put cov_type in kwargs so subclasses can decide in fit whether to 

1287 # use this generic implementation 

1288 if 'use_t' in kwargs: 

1289 use_t = kwargs['use_t'] 

1290 self.use_t = use_t if use_t is not None else False 

1291 if 'cov_type' in kwargs: 

1292 cov_type = kwargs.get('cov_type', 'nonrobust') 

1293 cov_kwds = kwargs.get('cov_kwds', {}) 

1294 

1295 if cov_type == 'nonrobust': 

1296 self.cov_type = 'nonrobust' 

1297 self.cov_kwds = {'description': 'Standard Errors assume that the ' + 

1298 'covariance matrix of the errors is correctly ' + 

1299 'specified.'} 

1300 else: 

1301 from statsmodels.base.covtype import get_robustcov_results 

1302 if cov_kwds is None: 

1303 cov_kwds = {} 

1304 use_t = self.use_t 

1305 # TODO: we should not need use_t in get_robustcov_results 

1306 get_robustcov_results(self, cov_type=cov_type, use_self=True, 

1307 use_t=use_t, **cov_kwds) 

1308 

1309 def normalized_cov_params(self): 

1310 """See specific model class docstring""" 

1311 raise NotImplementedError 

1312 

1313 def _get_robustcov_results(self, cov_type='nonrobust', use_self=True, 

1314 use_t=None, **cov_kwds): 

1315 if use_self is False: 

1316 raise ValueError("use_self should have been removed long ago. " 

1317 "See GH#4401") 

1318 from statsmodels.base.covtype import get_robustcov_results 

1319 if cov_kwds is None: 

1320 cov_kwds = {} 

1321 

1322 if cov_type == 'nonrobust': 

1323 self.cov_type = 'nonrobust' 

1324 self.cov_kwds = {'description': 'Standard Errors assume that the ' + 

1325 'covariance matrix of the errors is correctly ' + 

1326 'specified.'} 

1327 else: 

1328 # TODO: we should not need use_t in get_robustcov_results 

1329 get_robustcov_results(self, cov_type=cov_type, use_self=True, 

1330 use_t=use_t, **cov_kwds) 

1331 @property 

1332 def use_t(self): 

1333 """Flag indicating to use the Student's distribution in inference.""" 

1334 return self._use_t 

1335 

1336 @use_t.setter 

1337 def use_t(self, value): 

1338 self._use_t = bool(value) 

1339 

1340 @cached_value 

1341 def llf(self): 

1342 """Log-likelihood of model""" 

1343 return self.model.loglike(self.params) 

1344 

1345 @cached_value 

1346 def bse(self): 

1347 """The standard errors of the parameter estimates.""" 

1348 # Issue 3299 

1349 if ((not hasattr(self, 'cov_params_default')) and 

1350 (self.normalized_cov_params is None)): 

1351 bse_ = np.empty(len(self.params)) 

1352 bse_[:] = np.nan 

1353 else: 

1354 bse_ = np.sqrt(np.diag(self.cov_params())) 

1355 return bse_ 

1356 

1357 @cached_value 

1358 def tvalues(self): 

1359 """ 

1360 Return the t-statistic for a given parameter estimate. 

1361 """ 

1362 return self.params / self.bse 

1363 

1364 @cached_value 

1365 def pvalues(self): 

1366 """The two-tailed p values for the t-stats of the params.""" 

1367 if self.use_t: 

1368 df_resid = getattr(self, 'df_resid_inference', self.df_resid) 

1369 return stats.t.sf(np.abs(self.tvalues), df_resid) * 2 

1370 else: 

1371 return stats.norm.sf(np.abs(self.tvalues)) * 2 

1372 

1373 def cov_params(self, r_matrix=None, column=None, scale=None, cov_p=None, 

1374 other=None): 

1375 """ 

1376 Compute the variance/covariance matrix. 

1377 

1378 The variance/covariance matrix can be of a linear contrast of the 

1379 estimated parameters or all params multiplied by scale which will 

1380 usually be an estimate of sigma^2. Scale is assumed to be a scalar. 

1381 

1382 Parameters 

1383 ---------- 

1384 r_matrix : array_like 

1385 Can be 1d, or 2d. Can be used alone or with other. 

1386 column : array_like, optional 

1387 Must be used on its own. Can be 0d or 1d see below. 

1388 scale : float, optional 

1389 Can be specified or not. Default is None, which means that 

1390 the scale argument is taken from the model. 

1391 cov_p : ndarray, optional 

1392 The covariance of the parameters. If not provided, this value is 

1393 read from `self.normalized_cov_params` or 

1394 `self.cov_params_default`. 

1395 other : array_like, optional 

1396 Can be used when r_matrix is specified. 

1397 

1398 Returns 

1399 ------- 

1400 ndarray 

1401 The covariance matrix of the parameter estimates or of linear 

1402 combination of parameter estimates. See Notes. 

1403 

1404 Notes 

1405 ----- 

1406 (The below are assumed to be in matrix notation.) 

1407 

1408 If no argument is specified returns the covariance matrix of a model 

1409 ``(scale)*(X.T X)^(-1)`` 

1410 

1411 If contrast is specified it pre and post-multiplies as follows 

1412 ``(scale) * r_matrix (X.T X)^(-1) r_matrix.T`` 

1413 

1414 If contrast and other are specified returns 

1415 ``(scale) * r_matrix (X.T X)^(-1) other.T`` 

1416 

1417 If column is specified returns 

1418 ``(scale) * (X.T X)^(-1)[column,column]`` if column is 0d 

1419 

1420 OR 

1421 

1422 ``(scale) * (X.T X)^(-1)[column][:,column]`` if column is 1d 

1423 """ 

1424 if (hasattr(self, 'mle_settings') and 

1425 self.mle_settings['optimizer'] in ['l1', 'l1_cvxopt_cp']): 

1426 dot_fun = nan_dot 

1427 else: 

1428 dot_fun = np.dot 

1429 

1430 if (cov_p is None and self.normalized_cov_params is None and 

1431 not hasattr(self, 'cov_params_default')): 

1432 raise ValueError('need covariance of parameters for computing ' 

1433 '(unnormalized) covariances') 

1434 if column is not None and (r_matrix is not None or other is not None): 

1435 raise ValueError('Column should be specified without other ' 

1436 'arguments.') 

1437 if other is not None and r_matrix is None: 

1438 raise ValueError('other can only be specified with r_matrix') 

1439 

1440 if cov_p is None: 

1441 if hasattr(self, 'cov_params_default'): 

1442 cov_p = self.cov_params_default 

1443 else: 

1444 if scale is None: 

1445 scale = self.scale 

1446 cov_p = self.normalized_cov_params * scale 

1447 

1448 if column is not None: 

1449 column = np.asarray(column) 

1450 if column.shape == (): 

1451 return cov_p[column, column] 

1452 else: 

1453 return cov_p[column[:, None], column] 

1454 elif r_matrix is not None: 

1455 r_matrix = np.asarray(r_matrix) 

1456 if r_matrix.shape == (): 

1457 raise ValueError("r_matrix should be 1d or 2d") 

1458 if other is None: 

1459 other = r_matrix 

1460 else: 

1461 other = np.asarray(other) 

1462 tmp = dot_fun(r_matrix, dot_fun(cov_p, np.transpose(other))) 

1463 return tmp 

1464 else: # if r_matrix is None and column is None: 

1465 return cov_p 

1466 

1467 # TODO: make sure this works as needed for GLMs 

1468 def t_test(self, r_matrix, cov_p=None, scale=None, use_t=None): 

1469 """ 

1470 Compute a t-test for a each linear hypothesis of the form Rb = q. 

1471 

1472 Parameters 

1473 ---------- 

1474 r_matrix : {array_like, str, tuple} 

1475 One of: 

1476 

1477 - array : If an array is given, a p x k 2d array or length k 1d 

1478 array specifying the linear restrictions. It is assumed 

1479 that the linear combination is equal to zero. 

1480 - str : The full hypotheses to test can be given as a string. 

1481 See the examples. 

1482 - tuple : A tuple of arrays in the form (R, q). If q is given, 

1483 can be either a scalar or a length p row vector. 

1484 

1485 cov_p : array_like, optional 

1486 An alternative estimate for the parameter covariance matrix. 

1487 If None is given, self.normalized_cov_params is used. 

1488 scale : float, optional 

1489 An optional `scale` to use. Default is the scale specified 

1490 by the model fit. 

1491 

1492 .. deprecated:: 0.10.0 

1493 

1494 use_t : bool, optional 

1495 If use_t is None, then the default of the model is used. If use_t 

1496 is True, then the p-values are based on the t distribution. If 

1497 use_t is False, then the p-values are based on the normal 

1498 distribution. 

1499 

1500 Returns 

1501 ------- 

1502 ContrastResults 

1503 The results for the test are attributes of this results instance. 

1504 The available results have the same elements as the parameter table 

1505 in `summary()`. 

1506 

1507 See Also 

1508 -------- 

1509 tvalues : Individual t statistics for the estimated parameters. 

1510 f_test : Perform an F tests on model parameters. 

1511 patsy.DesignInfo.linear_constraint : Specify a linear constraint. 

1512 

1513 Examples 

1514 -------- 

1515 >>> import numpy as np 

1516 >>> import statsmodels.api as sm 

1517 >>> data = sm.datasets.longley.load(as_pandas=False) 

1518 >>> data.exog = sm.add_constant(data.exog) 

1519 >>> results = sm.OLS(data.endog, data.exog).fit() 

1520 >>> r = np.zeros_like(results.params) 

1521 >>> r[5:] = [1,-1] 

1522 >>> print(r) 

1523 [ 0. 0. 0. 0. 0. 1. -1.] 

1524 

1525 r tests that the coefficients on the 5th and 6th independent 

1526 variable are the same. 

1527 

1528 >>> T_test = results.t_test(r) 

1529 >>> print(T_test) 

1530 Test for Constraints 

1531 ============================================================================== 

1532 coef std err t P>|t| [0.025 0.975] 

1533 ------------------------------------------------------------------------------ 

1534 c0 -1829.2026 455.391 -4.017 0.003 -2859.368 -799.037 

1535 ============================================================================== 

1536 >>> T_test.effect 

1537 -1829.2025687192481 

1538 >>> T_test.sd 

1539 455.39079425193762 

1540 >>> T_test.tvalue 

1541 -4.0167754636411717 

1542 >>> T_test.pvalue 

1543 0.0015163772380899498 

1544 

1545 Alternatively, you can specify the hypothesis tests using a string 

1546 

1547 >>> from statsmodels.formula.api import ols 

1548 >>> dta = sm.datasets.longley.load_pandas().data 

1549 >>> formula = 'TOTEMP ~ GNPDEFL + GNP + UNEMP + ARMED + POP + YEAR' 

1550 >>> results = ols(formula, dta).fit() 

1551 >>> hypotheses = 'GNPDEFL = GNP, UNEMP = 2, YEAR/1829 = 1' 

1552 >>> t_test = results.t_test(hypotheses) 

1553 >>> print(t_test) 

1554 Test for Constraints 

1555 ============================================================================== 

1556 coef std err t P>|t| [0.025 0.975] 

1557 ------------------------------------------------------------------------------ 

1558 c0 15.0977 84.937 0.178 0.863 -177.042 207.238 

1559 c1 -2.0202 0.488 -8.231 0.000 -3.125 -0.915 

1560 c2 1.0001 0.249 0.000 1.000 0.437 1.563 

1561 ============================================================================== 

1562 """ 

1563 if scale is not None: 

1564 import warnings 

1565 warnings.warn('scale is has no effect and is deprecated. It will' 

1566 'be removed in the next version.', 

1567 DeprecationWarning) 

1568 

1569 from patsy import DesignInfo 

1570 names = self.model.data.cov_names 

1571 LC = DesignInfo(names).linear_constraint(r_matrix) 

1572 r_matrix, q_matrix = LC.coefs, LC.constants 

1573 num_ttests = r_matrix.shape[0] 

1574 num_params = r_matrix.shape[1] 

1575 

1576 if (cov_p is None and self.normalized_cov_params is None and 

1577 not hasattr(self, 'cov_params_default')): 

1578 raise ValueError('Need covariance of parameters for computing ' 

1579 'T statistics') 

1580 params = self.params.ravel() 

1581 if num_params != params.shape[0]: 

1582 raise ValueError('r_matrix and params are not aligned') 

1583 if q_matrix is None: 

1584 q_matrix = np.zeros(num_ttests) 

1585 else: 

1586 q_matrix = np.asarray(q_matrix) 

1587 q_matrix = q_matrix.squeeze() 

1588 if q_matrix.size > 1: 

1589 if q_matrix.shape[0] != num_ttests: 

1590 raise ValueError("r_matrix and q_matrix must have the same " 

1591 "number of rows") 

1592 

1593 if use_t is None: 

1594 # switch to use_t false if undefined 

1595 use_t = (hasattr(self, 'use_t') and self.use_t) 

1596 

1597 _effect = np.dot(r_matrix, params) 

1598 

1599 # Perform the test 

1600 if num_ttests > 1: 

1601 _sd = np.sqrt(np.diag(self.cov_params( 

1602 r_matrix=r_matrix, cov_p=cov_p))) 

1603 else: 

1604 _sd = np.sqrt(self.cov_params(r_matrix=r_matrix, cov_p=cov_p)) 

1605 _t = (_effect - q_matrix) * recipr(_sd) 

1606 

1607 df_resid = getattr(self, 'df_resid_inference', self.df_resid) 

1608 

1609 if use_t: 

1610 return ContrastResults(effect=_effect, t=_t, sd=_sd, 

1611 df_denom=df_resid) 

1612 else: 

1613 return ContrastResults(effect=_effect, statistic=_t, sd=_sd, 

1614 df_denom=df_resid, 

1615 distribution='norm') 

1616 

1617 def f_test(self, r_matrix, cov_p=None, scale=1.0, invcov=None): 

1618 """ 

1619 Compute the F-test for a joint linear hypothesis. 

1620 

1621 This is a special case of `wald_test` that always uses the F 

1622 distribution. 

1623 

1624 Parameters 

1625 ---------- 

1626 r_matrix : {array_like, str, tuple} 

1627 One of: 

1628 

1629 - array : An r x k array where r is the number of restrictions to 

1630 test and k is the number of regressors. It is assumed 

1631 that the linear combination is equal to zero. 

1632 - str : The full hypotheses to test can be given as a string. 

1633 See the examples. 

1634 - tuple : A tuple of arrays in the form (R, q), ``q`` can be 

1635 either a scalar or a length k row vector. 

1636 

1637 cov_p : array_like, optional 

1638 An alternative estimate for the parameter covariance matrix. 

1639 If None is given, self.normalized_cov_params is used. 

1640 scale : float, optional 

1641 Default is 1.0 for no scaling. 

1642 

1643 .. deprecated:: 0.10.0 

1644 

1645 invcov : array_like, optional 

1646 A q x q array to specify an inverse covariance matrix based on a 

1647 restrictions matrix. 

1648 

1649 Returns 

1650 ------- 

1651 ContrastResults 

1652 The results for the test are attributes of this results instance. 

1653 

1654 See Also 

1655 -------- 

1656 t_test : Perform a single hypothesis test. 

1657 wald_test : Perform a Wald-test using a quadratic form. 

1658 statsmodels.stats.contrast.ContrastResults : Test results. 

1659 patsy.DesignInfo.linear_constraint : Specify a linear constraint. 

1660 

1661 Notes 

1662 ----- 

1663 The matrix `r_matrix` is assumed to be non-singular. More precisely, 

1664 

1665 r_matrix (pX pX.T) r_matrix.T 

1666 

1667 is assumed invertible. Here, pX is the generalized inverse of the 

1668 design matrix of the model. There can be problems in non-OLS models 

1669 where the rank of the covariance of the noise is not full. 

1670 

1671 Examples 

1672 -------- 

1673 >>> import numpy as np 

1674 >>> import statsmodels.api as sm 

1675 >>> data = sm.datasets.longley.load(as_pandas=False) 

1676 >>> data.exog = sm.add_constant(data.exog) 

1677 >>> results = sm.OLS(data.endog, data.exog).fit() 

1678 >>> A = np.identity(len(results.params)) 

1679 >>> A = A[1:,:] 

1680 

1681 This tests that each coefficient is jointly statistically 

1682 significantly different from zero. 

1683 

1684 >>> print(results.f_test(A)) 

1685 <F test: F=array([[ 330.28533923]]), p=4.984030528700946e-10, df_denom=9, df_num=6> 

1686 

1687 Compare this to 

1688 

1689 >>> results.fvalue 

1690 330.2853392346658 

1691 >>> results.f_pvalue 

1692 4.98403096572e-10 

1693 

1694 >>> B = np.array(([0,0,1,-1,0,0,0],[0,0,0,0,0,1,-1])) 

1695 

1696 This tests that the coefficient on the 2nd and 3rd regressors are 

1697 equal and jointly that the coefficient on the 5th and 6th regressors 

1698 are equal. 

1699 

1700 >>> print(results.f_test(B)) 

1701 <F test: F=array([[ 9.74046187]]), p=0.005605288531708235, df_denom=9, df_num=2> 

1702 

1703 Alternatively, you can specify the hypothesis tests using a string 

1704 

1705 >>> from statsmodels.datasets import longley 

1706 >>> from statsmodels.formula.api import ols 

1707 >>> dta = longley.load_pandas().data 

1708 >>> formula = 'TOTEMP ~ GNPDEFL + GNP + UNEMP + ARMED + POP + YEAR' 

1709 >>> results = ols(formula, dta).fit() 

1710 >>> hypotheses = '(GNPDEFL = GNP), (UNEMP = 2), (YEAR/1829 = 1)' 

1711 >>> f_test = results.f_test(hypotheses) 

1712 >>> print(f_test) 

1713 <F test: F=array([[ 144.17976065]]), p=6.322026217355609e-08, df_denom=9, df_num=3> 

1714 """ 

1715 if scale != 1.0: 

1716 import warnings 

1717 warnings.warn('scale is has no effect and is deprecated. It will' 

1718 'be removed in the next version.', 

1719 DeprecationWarning) 

1720 

1721 res = self.wald_test(r_matrix, cov_p=cov_p, invcov=invcov, use_f=True) 

1722 return res 

1723 

1724 # TODO: untested for GLMs? 

1725 def wald_test(self, r_matrix, cov_p=None, scale=1.0, invcov=None, 

1726 use_f=None, df_constraints=None): 

1727 """ 

1728 Compute a Wald-test for a joint linear hypothesis. 

1729 

1730 Parameters 

1731 ---------- 

1732 r_matrix : {array_like, str, tuple} 

1733 One of: 

1734 

1735 - array : An r x k array where r is the number of restrictions to 

1736 test and k is the number of regressors. It is assumed that the 

1737 linear combination is equal to zero. 

1738 - str : The full hypotheses to test can be given as a string. 

1739 See the examples. 

1740 - tuple : A tuple of arrays in the form (R, q), ``q`` can be 

1741 either a scalar or a length p row vector. 

1742 

1743 cov_p : array_like, optional 

1744 An alternative estimate for the parameter covariance matrix. 

1745 If None is given, self.normalized_cov_params is used. 

1746 scale : float, optional 

1747 Default is 1.0 for no scaling. 

1748 

1749 .. deprecated:: 0.10.0 

1750 

1751 invcov : array_like, optional 

1752 A q x q array to specify an inverse covariance matrix based on a 

1753 restrictions matrix. 

1754 use_f : bool 

1755 If True, then the F-distribution is used. If False, then the 

1756 asymptotic distribution, chisquare is used. If use_f is None, then 

1757 the F distribution is used if the model specifies that use_t is True. 

1758 The test statistic is proportionally adjusted for the distribution 

1759 by the number of constraints in the hypothesis. 

1760 df_constraints : int, optional 

1761 The number of constraints. If not provided the number of 

1762 constraints is determined from r_matrix. 

1763 

1764 Returns 

1765 ------- 

1766 ContrastResults 

1767 The results for the test are attributes of this results instance. 

1768 

1769 See Also 

1770 -------- 

1771 f_test : Perform an F tests on model parameters. 

1772 t_test : Perform a single hypothesis test. 

1773 statsmodels.stats.contrast.ContrastResults : Test results. 

1774 patsy.DesignInfo.linear_constraint : Specify a linear constraint. 

1775 

1776 Notes 

1777 ----- 

1778 The matrix `r_matrix` is assumed to be non-singular. More precisely, 

1779 

1780 r_matrix (pX pX.T) r_matrix.T 

1781 

1782 is assumed invertible. Here, pX is the generalized inverse of the 

1783 design matrix of the model. There can be problems in non-OLS models 

1784 where the rank of the covariance of the noise is not full. 

1785 """ 

1786 if scale != 1.0: 

1787 import warnings 

1788 warnings.warn('scale is has no effect and is deprecated. It will' 

1789 'be removed in the next version.', 

1790 DeprecationWarning) 

1791 

1792 if use_f is None: 

1793 # switch to use_t false if undefined 

1794 use_f = (hasattr(self, 'use_t') and self.use_t) 

1795 

1796 from patsy import DesignInfo 

1797 names = self.model.data.cov_names 

1798 params = self.params.ravel() 

1799 LC = DesignInfo(names).linear_constraint(r_matrix) 

1800 r_matrix, q_matrix = LC.coefs, LC.constants 

1801 

1802 if (self.normalized_cov_params is None and cov_p is None and 

1803 invcov is None and not hasattr(self, 'cov_params_default')): 

1804 raise ValueError('need covariance of parameters for computing ' 

1805 'F statistics') 

1806 

1807 cparams = np.dot(r_matrix, params[:, None]) 

1808 J = float(r_matrix.shape[0]) # number of restrictions 

1809 

1810 if q_matrix is None: 

1811 q_matrix = np.zeros(J) 

1812 else: 

1813 q_matrix = np.asarray(q_matrix) 

1814 if q_matrix.ndim == 1: 

1815 q_matrix = q_matrix[:, None] 

1816 if q_matrix.shape[0] != J: 

1817 raise ValueError("r_matrix and q_matrix must have the same " 

1818 "number of rows") 

1819 Rbq = cparams - q_matrix 

1820 if invcov is None: 

1821 cov_p = self.cov_params(r_matrix=r_matrix, cov_p=cov_p) 

1822 if np.isnan(cov_p).max(): 

1823 raise ValueError("r_matrix performs f_test for using " 

1824 "dimensions that are asymptotically " 

1825 "non-normal") 

1826 invcov = np.linalg.pinv(cov_p) 

1827 J_ = np.linalg.matrix_rank(cov_p) 

1828 if J_ < J: 

1829 import warnings 

1830 warnings.warn('covariance of constraints does not have full ' 

1831 'rank. The number of constraints is %d, but ' 

1832 'rank is %d' % (J, J_), ValueWarning) 

1833 J = J_ 

1834 

1835 # TODO streamline computation, we do not need to compute J if given 

1836 if df_constraints is not None: 

1837 # let caller override J by df_constraint 

1838 J = df_constraints 

1839 

1840 if (hasattr(self, 'mle_settings') and 

1841 self.mle_settings['optimizer'] in ['l1', 'l1_cvxopt_cp']): 

1842 F = nan_dot(nan_dot(Rbq.T, invcov), Rbq) 

1843 else: 

1844 F = np.dot(np.dot(Rbq.T, invcov), Rbq) 

1845 

1846 df_resid = getattr(self, 'df_resid_inference', self.df_resid) 

1847 if use_f: 

1848 F /= J 

1849 return ContrastResults(F=F, df_denom=df_resid, 

1850 df_num=J) #invcov.shape[0]) 

1851 else: 

1852 return ContrastResults(chi2=F, df_denom=J, statistic=F, 

1853 distribution='chi2', distargs=(J,)) 

1854 

1855 def wald_test_terms(self, skip_single=False, extra_constraints=None, 

1856 combine_terms=None): 

1857 """ 

1858 Compute a sequence of Wald tests for terms over multiple columns. 

1859 

1860 This computes joined Wald tests for the hypothesis that all 

1861 coefficients corresponding to a `term` are zero. 

1862 `Terms` are defined by the underlying formula or by string matching. 

1863 

1864 Parameters 

1865 ---------- 

1866 skip_single : bool 

1867 If true, then terms that consist only of a single column and, 

1868 therefore, refers only to a single parameter is skipped. 

1869 If false, then all terms are included. 

1870 extra_constraints : ndarray 

1871 Additional constraints to test. Note that this input has not been 

1872 tested. 

1873 combine_terms : {list[str], None} 

1874 Each string in this list is matched to the name of the terms or 

1875 the name of the exogenous variables. All columns whose name 

1876 includes that string are combined in one joint test. 

1877 

1878 Returns 

1879 ------- 

1880 WaldTestResults 

1881 The result instance contains `table` which is a pandas DataFrame 

1882 with the test results: test statistic, degrees of freedom and 

1883 pvalues. 

1884 

1885 Examples 

1886 -------- 

1887 >>> res_ols = ols("np.log(Days+1) ~ C(Duration, Sum)*C(Weight, Sum)", data).fit() 

1888 >>> res_ols.wald_test_terms() 

1889 <class 'statsmodels.stats.contrast.WaldTestResults'> 

1890 F P>F df constraint df denom 

1891 Intercept 279.754525 2.37985521351e-22 1 51 

1892 C(Duration, Sum) 5.367071 0.0245738436636 1 51 

1893 C(Weight, Sum) 12.432445 3.99943118767e-05 2 51 

1894 C(Duration, Sum):C(Weight, Sum) 0.176002 0.83912310946 2 51 

1895 

1896 >>> res_poi = Poisson.from_formula("Days ~ C(Weight) * C(Duration)", \ 

1897 data).fit(cov_type='HC0') 

1898 >>> wt = res_poi.wald_test_terms(skip_single=False, \ 

1899 combine_terms=['Duration', 'Weight']) 

1900 >>> print(wt) 

1901 chi2 P>chi2 df constraint 

1902 Intercept 15.695625 7.43960374424e-05 1 

1903 C(Weight) 16.132616 0.000313940174705 2 

1904 C(Duration) 1.009147 0.315107378931 1 

1905 C(Weight):C(Duration) 0.216694 0.897315972824 2 

1906 Duration 11.187849 0.010752286833 3 

1907 Weight 30.263368 4.32586407145e-06 4 

1908 """ 

1909 # lazy import 

1910 from collections import defaultdict 

1911 

1912 result = self 

1913 if extra_constraints is None: 

1914 extra_constraints = [] 

1915 if combine_terms is None: 

1916 combine_terms = [] 

1917 design_info = getattr(result.model.data, 'design_info', None) 

1918 

1919 if design_info is None and extra_constraints is None: 

1920 raise ValueError('no constraints, nothing to do') 

1921 

1922 identity = np.eye(len(result.params)) 

1923 constraints = [] 

1924 combined = defaultdict(list) 

1925 if design_info is not None: 

1926 for term in design_info.terms: 

1927 cols = design_info.slice(term) 

1928 name = term.name() 

1929 constraint_matrix = identity[cols] 

1930 

1931 # check if in combined 

1932 for cname in combine_terms: 

1933 if cname in name: 

1934 combined[cname].append(constraint_matrix) 

1935 

1936 k_constraint = constraint_matrix.shape[0] 

1937 if skip_single: 

1938 if k_constraint == 1: 

1939 continue 

1940 

1941 constraints.append((name, constraint_matrix)) 

1942 

1943 combined_constraints = [] 

1944 for cname in combine_terms: 

1945 combined_constraints.append((cname, np.vstack(combined[cname]))) 

1946 else: 

1947 # check by exog/params names if there is no formula info 

1948 for col, name in enumerate(result.model.exog_names): 

1949 constraint_matrix = np.atleast_2d(identity[col]) 

1950 

1951 # check if in combined 

1952 for cname in combine_terms: 

1953 if cname in name: 

1954 combined[cname].append(constraint_matrix) 

1955 

1956 if skip_single: 

1957 continue 

1958 

1959 constraints.append((name, constraint_matrix)) 

1960 

1961 combined_constraints = [] 

1962 for cname in combine_terms: 

1963 combined_constraints.append((cname, np.vstack(combined[cname]))) 

1964 

1965 use_t = result.use_t 

1966 distribution = ['chi2', 'F'][use_t] 

1967 

1968 res_wald = [] 

1969 index = [] 

1970 for name, constraint in constraints + combined_constraints + extra_constraints: 

1971 wt = result.wald_test(constraint) 

1972 row = [wt.statistic.item(), wt.pvalue.item(), constraint.shape[0]] 

1973 if use_t: 

1974 row.append(wt.df_denom) 

1975 res_wald.append(row) 

1976 index.append(name) 

1977 

1978 # distribution nerutral names 

1979 col_names = ['statistic', 'pvalue', 'df_constraint'] 

1980 if use_t: 

1981 col_names.append('df_denom') 

1982 # TODO: maybe move DataFrame creation to results class 

1983 from pandas import DataFrame 

1984 table = DataFrame(res_wald, index=index, columns=col_names) 

1985 res = WaldTestResults(None, distribution, None, table=table) 

1986 # TODO: remove temp again, added for testing 

1987 res.temp = constraints + combined_constraints + extra_constraints 

1988 return res 

1989 

1990 def t_test_pairwise(self, term_name, method='hs', alpha=0.05, 

1991 factor_labels=None): 

1992 """ 

1993 Perform pairwise t_test with multiple testing corrected p-values. 

1994 

1995 This uses the formula design_info encoding contrast matrix and should 

1996 work for all encodings of a main effect. 

1997 

1998 Parameters 

1999 ---------- 

2000 term_name : str 

2001 The name of the term for which pairwise comparisons are computed. 

2002 Term names for categorical effects are created by patsy and 

2003 correspond to the main part of the exog names. 

2004 method : {str, list[str]} 

2005 The multiple testing p-value correction to apply. The default is 

2006 'hs'. See stats.multipletesting. 

2007 alpha : float 

2008 The significance level for multiple testing reject decision. 

2009 factor_labels : {list[str], None} 

2010 Labels for the factor levels used for pairwise labels. If not 

2011 provided, then the labels from the formula design_info are used. 

2012 

2013 Returns 

2014 ------- 

2015 MultiCompResult 

2016 The results are stored as attributes, the main attributes are the 

2017 following two. Other attributes are added for debugging purposes 

2018 or as background information. 

2019 

2020 - result_frame : pandas DataFrame with t_test results and multiple 

2021 testing corrected p-values. 

2022 - contrasts : matrix of constraints of the null hypothesis in the 

2023 t_test. 

2024 

2025 Notes 

2026 ----- 

2027 Status: experimental. Currently only checked for treatment coding with 

2028 and without specified reference level. 

2029 

2030 Currently there are no multiple testing corrected confidence intervals 

2031 available. 

2032 

2033 Examples 

2034 -------- 

2035 >>> res = ols("np.log(Days+1) ~ C(Weight) + C(Duration)", data).fit() 

2036 >>> pw = res.t_test_pairwise("C(Weight)") 

2037 >>> pw.result_frame 

2038 coef std err t P>|t| Conf. Int. Low 

2039 2-1 0.632315 0.230003 2.749157 8.028083e-03 0.171563 

2040 3-1 1.302555 0.230003 5.663201 5.331513e-07 0.841803 

2041 3-2 0.670240 0.230003 2.914044 5.119126e-03 0.209488 

2042 Conf. Int. Upp. pvalue-hs reject-hs 

2043 2-1 1.093067 0.010212 True 

2044 3-1 1.763307 0.000002 True 

2045 3-2 1.130992 0.010212 True 

2046 """ 

2047 res = t_test_pairwise(self, term_name, method=method, alpha=alpha, 

2048 factor_labels=factor_labels) 

2049 return res 

2050 

2051 def conf_int(self, alpha=.05, cols=None): 

2052 """ 

2053 Construct confidence interval for the fitted parameters. 

2054 

2055 Parameters 

2056 ---------- 

2057 alpha : float, optional 

2058 The significance level for the confidence interval. The default 

2059 `alpha` = .05 returns a 95% confidence interval. 

2060 cols : array_like, optional 

2061 Specifies which confidence intervals to return. 

2062 

2063 Returns 

2064 ------- 

2065 array_like 

2066 Each row contains [lower, upper] limits of the confidence interval 

2067 for the corresponding parameter. The first column contains all 

2068 lower, the second column contains all upper limits. 

2069 

2070 Notes 

2071 ----- 

2072 The confidence interval is based on the standard normal distribution 

2073 if self.use_t is False. If self.use_t is True, then uses a Student's t 

2074 with self.df_resid_inference (or self.df_resid if df_resid_inference is 

2075 not defined) degrees of freedom. 

2076 

2077 Examples 

2078 -------- 

2079 >>> import statsmodels.api as sm 

2080 >>> data = sm.datasets.longley.load(as_pandas=False) 

2081 >>> data.exog = sm.add_constant(data.exog) 

2082 >>> results = sm.OLS(data.endog, data.exog).fit() 

2083 >>> results.conf_int() 

2084 array([[-5496529.48322745, -1467987.78596704], 

2085 [ -177.02903529, 207.15277984], 

2086 [ -0.1115811 , 0.03994274], 

2087 [ -3.12506664, -0.91539297], 

2088 [ -1.5179487 , -0.54850503], 

2089 [ -0.56251721, 0.460309 ], 

2090 [ 798.7875153 , 2859.51541392]]) 

2091 

2092 >>> results.conf_int(cols=(2,3)) 

2093 array([[-0.1115811 , 0.03994274], 

2094 [-3.12506664, -0.91539297]]) 

2095 """ 

2096 bse = self.bse 

2097 

2098 if self.use_t: 

2099 dist = stats.t 

2100 df_resid = getattr(self, 'df_resid_inference', self.df_resid) 

2101 q = dist.ppf(1 - alpha / 2, df_resid) 

2102 else: 

2103 dist = stats.norm 

2104 q = dist.ppf(1 - alpha / 2) 

2105 

2106 params = self.params 

2107 lower = params - q * bse 

2108 upper = params + q * bse 

2109 if cols is not None: 

2110 cols = np.asarray(cols) 

2111 lower = lower[cols] 

2112 upper = upper[cols] 

2113 return np.asarray(lzip(lower, upper)) 

2114 

2115 def save(self, fname, remove_data=False): 

2116 """ 

2117 Save a pickle of this instance. 

2118 

2119 Parameters 

2120 ---------- 

2121 fname : {str, handle} 

2122 A string filename or a file handle. 

2123 remove_data : bool 

2124 If False (default), then the instance is pickled without changes. 

2125 If True, then all arrays with length nobs are set to None before 

2126 pickling. See the remove_data method. 

2127 In some cases not all arrays will be set to None. 

2128 

2129 Notes 

2130 ----- 

2131 If remove_data is true and the model result does not implement a 

2132 remove_data method then this will raise an exception. 

2133 """ 

2134 

2135 from statsmodels.iolib.smpickle import save_pickle 

2136 

2137 if remove_data: 

2138 self.remove_data() 

2139 

2140 save_pickle(self, fname) 

2141 

2142 @classmethod 

2143 def load(cls, fname): 

2144 """ 

2145 Load a pickled results instance 

2146 

2147 .. warning:: 

2148 

2149 Loading pickled models is not secure against erroneous or 

2150 maliciously constructed data. Never unpickle data received from 

2151 an untrusted or unauthenticated source. 

2152 

2153 Parameters 

2154 ---------- 

2155 fname : {str, handle} 

2156 A string filename or a file handle. 

2157 

2158 Returns 

2159 ------- 

2160 Results 

2161 The unpickled results instance. 

2162 """ 

2163 

2164 from statsmodels.iolib.smpickle import load_pickle 

2165 return load_pickle(fname) 

2166 

2167 def remove_data(self): 

2168 """ 

2169 Remove data arrays, all nobs arrays from result and model. 

2170 

2171 This reduces the size of the instance, so it can be pickled with less 

2172 memory. Currently tested for use with predict from an unpickled 

2173 results and model instance. 

2174 

2175 .. warning:: 

2176 

2177 Since data and some intermediate results have been removed 

2178 calculating new statistics that require them will raise exceptions. 

2179 The exception will occur the first time an attribute is accessed 

2180 that has been set to None. 

2181 

2182 Not fully tested for time series models, tsa, and might delete too much 

2183 for prediction or not all that would be possible. 

2184 

2185 The lists of arrays to delete are maintained as attributes of 

2186 the result and model instance, except for cached values. These 

2187 lists could be changed before calling remove_data. 

2188 

2189 The attributes to remove are named in: 

2190 

2191 model._data_attr : arrays attached to both the model instance 

2192 and the results instance with the same attribute name. 

2193 

2194 result.data_in_cache : arrays that may exist as values in 

2195 result._cache (TODO : should privatize name) 

2196 

2197 result._data_attr_model : arrays attached to the model 

2198 instance but not to the results instance 

2199 """ 

2200 cls = self.__class__ 

2201 # Note: we cannot just use `getattr(cls, x)` or `getattr(self, x)` 

2202 # because of redirection involved with property-like accessors 

2203 cls_attrs = {} 

2204 for name in dir(cls): 

2205 try: 

2206 attr = object.__getattribute__(cls, name) 

2207 except AttributeError: 

2208 pass 

2209 else: 

2210 cls_attrs[name] = attr 

2211 data_attrs = [x for x in cls_attrs 

2212 if isinstance(cls_attrs[x], cached_data)] 

2213 value_attrs = [x for x in cls_attrs 

2214 if isinstance(cls_attrs[x], cached_value)] 

2215 # make sure the cached for value_attrs are evaluated; this needs to 

2216 # occur _before_ any other attributes are removed. 

2217 for name in value_attrs: 

2218 getattr(self, name) 

2219 for name in data_attrs: 

2220 self._cache[name] = None 

2221 

2222 def wipe(obj, att): 

2223 # get to last element in attribute path 

2224 p = att.split('.') 

2225 att_ = p.pop(-1) 

2226 try: 

2227 obj_ = reduce(getattr, [obj] + p) 

2228 if hasattr(obj_, att_): 

2229 setattr(obj_, att_, None) 

2230 except AttributeError: 

2231 pass 

2232 

2233 model_only = ['model.' + i for i in getattr(self, "_data_attr_model", [])] 

2234 model_attr = ['model.' + i for i in self.model._data_attr] 

2235 for att in self._data_attr + model_attr + model_only: 

2236 if att in data_attrs: 

2237 # these have been handled above, and trying to call wipe 

2238 # would raise an Exception anyway, so skip these 

2239 continue 

2240 wipe(self, att) 

2241 

2242 data_in_cache = getattr(self, 'data_in_cache', []) 

2243 data_in_cache += ['fittedvalues', 'resid', 'wresid'] 

2244 for key in data_in_cache: 

2245 try: 

2246 self._cache[key] = None 

2247 except (AttributeError, KeyError): 

2248 pass 

2249 

2250 

2251class LikelihoodResultsWrapper(wrap.ResultsWrapper): 

2252 _attrs = { 

2253 'params': 'columns', 

2254 'bse': 'columns', 

2255 'pvalues': 'columns', 

2256 'tvalues': 'columns', 

2257 'resid': 'rows', 

2258 'fittedvalues': 'rows', 

2259 'normalized_cov_params': 'cov', 

2260 } 

2261 

2262 _wrap_attrs = _attrs 

2263 _wrap_methods = { 

2264 'cov_params': 'cov', 

2265 'conf_int': 'columns' 

2266 } 

2267 

2268wrap.populate_wrapper(LikelihoodResultsWrapper, # noqa:E305 

2269 LikelihoodModelResults) 

2270 

2271 

2272class ResultMixin(object): 

2273 

2274 @cache_readonly 

2275 def df_modelwc(self): 

2276 """Model WC""" 

2277 # collect different ways of defining the number of parameters, used for 

2278 # aic, bic 

2279 if hasattr(self, 'df_model'): 

2280 if hasattr(self, 'hasconst'): 

2281 hasconst = self.hasconst 

2282 else: 

2283 # default assumption 

2284 hasconst = 1 

2285 return self.df_model + hasconst 

2286 else: 

2287 return self.params.size 

2288 

2289 @cache_readonly 

2290 def aic(self): 

2291 """Akaike information criterion""" 

2292 return -2 * self.llf + 2 * (self.df_modelwc) 

2293 

2294 @cache_readonly 

2295 def bic(self): 

2296 """Bayesian information criterion""" 

2297 return -2 * self.llf + np.log(self.nobs) * (self.df_modelwc) 

2298 

2299 @cache_readonly 

2300 def score_obsv(self): 

2301 """cached Jacobian of log-likelihood 

2302 """ 

2303 return self.model.score_obs(self.params) 

2304 

2305 @cache_readonly 

2306 def hessv(self): 

2307 """cached Hessian of log-likelihood 

2308 """ 

2309 return self.model.hessian(self.params) 

2310 

2311 @cache_readonly 

2312 def covjac(self): 

2313 """ 

2314 covariance of parameters based on outer product of jacobian of 

2315 log-likelihood 

2316 """ 

2317 # if not hasattr(self, '_results'): 

2318 # raise ValueError('need to call fit first') 

2319 # #self.fit() 

2320 # self.jacv = jacv = self.jac(self._results.params) 

2321 jacv = self.score_obsv 

2322 return np.linalg.inv(np.dot(jacv.T, jacv)) 

2323 

2324 @cache_readonly 

2325 def covjhj(self): 

2326 """covariance of parameters based on HJJH 

2327 

2328 dot product of Hessian, Jacobian, Jacobian, Hessian of likelihood 

2329 

2330 name should be covhjh 

2331 """ 

2332 jacv = self.score_obsv 

2333 hessv = self.hessv 

2334 hessinv = np.linalg.inv(hessv) 

2335 # self.hessinv = hessin = self.cov_params() 

2336 return np.dot(hessinv, np.dot(np.dot(jacv.T, jacv), hessinv)) 

2337 

2338 @cache_readonly 

2339 def bsejhj(self): 

2340 """standard deviation of parameter estimates based on covHJH 

2341 """ 

2342 return np.sqrt(np.diag(self.covjhj)) 

2343 

2344 @cache_readonly 

2345 def bsejac(self): 

2346 """standard deviation of parameter estimates based on covjac 

2347 """ 

2348 return np.sqrt(np.diag(self.covjac)) 

2349 

2350 def bootstrap(self, nrep=100, method='nm', disp=0, store=1): 

2351 """simple bootstrap to get mean and variance of estimator 

2352 

2353 see notes 

2354 

2355 Parameters 

2356 ---------- 

2357 nrep : int 

2358 number of bootstrap replications 

2359 method : str 

2360 optimization method to use 

2361 disp : bool 

2362 If true, then optimization prints results 

2363 store : bool 

2364 If true, then parameter estimates for all bootstrap iterations 

2365 are attached in self.bootstrap_results 

2366 

2367 Returns 

2368 ------- 

2369 mean : ndarray 

2370 mean of parameter estimates over bootstrap replications 

2371 std : ndarray 

2372 standard deviation of parameter estimates over bootstrap 

2373 replications 

2374 

2375 Notes 

2376 ----- 

2377 This was mainly written to compare estimators of the standard errors of 

2378 the parameter estimates. It uses independent random sampling from the 

2379 original endog and exog, and therefore is only correct if observations 

2380 are independently distributed. 

2381 

2382 This will be moved to apply only to models with independently 

2383 distributed observations. 

2384 """ 

2385 results = [] 

2386 print(self.model.__class__) 

2387 hascloneattr = True if hasattr(self.model, 'cloneattr') else False 

2388 for i in range(nrep): 

2389 rvsind = np.random.randint(self.nobs, size=self.nobs) 

2390 # this needs to set startparam and get other defining attributes 

2391 # need a clone method on model 

2392 if self.exog is not None: 

2393 exog_resamp = self.exog[rvsind, :] 

2394 else: 

2395 exog_resamp = None 

2396 fitmod = self.model.__class__(self.endog[rvsind], 

2397 exog=exog_resamp) 

2398 if hascloneattr: 

2399 for attr in self.model.cloneattr: 

2400 setattr(fitmod, attr, getattr(self.model, attr)) 

2401 

2402 fitres = fitmod.fit(method=method, disp=disp) 

2403 results.append(fitres.params) 

2404 results = np.array(results) 

2405 if store: 

2406 self.bootstrap_results = results 

2407 return results.mean(0), results.std(0), results 

2408 

2409 def get_nlfun(self, fun): 

2410 """ 

2411 get_nlfun 

2412 

2413 This is not Implemented 

2414 """ 

2415 # I think this is supposed to get the delta method that is currently 

2416 # in miscmodels count (as part of Poisson example) 

2417 raise NotImplementedError 

2418 

2419 

2420class GenericLikelihoodModelResults(LikelihoodModelResults, ResultMixin): 

2421 """ 

2422 A results class for the discrete dependent variable models. 

2423 

2424 ..Warning : 

2425 

2426 The following description has not been updated to this version/class. 

2427 Where are AIC, BIC, ....? docstring looks like copy from discretemod 

2428 

2429 Parameters 

2430 ---------- 

2431 model : A DiscreteModel instance 

2432 mlefit : instance of LikelihoodResults 

2433 This contains the numerical optimization results as returned by 

2434 LikelihoodModel.fit(), in a superclass of GnericLikelihoodModels 

2435 

2436 

2437 Attributes 

2438 ---------- 

2439 aic : float 

2440 Akaike information criterion. -2*(`llf` - p) where p is the number 

2441 of regressors including the intercept. 

2442 bic : float 

2443 Bayesian information criterion. -2*`llf` + ln(`nobs`)*p where p is the 

2444 number of regressors including the intercept. 

2445 bse : ndarray 

2446 The standard errors of the coefficients. 

2447 df_resid : float 

2448 See model definition. 

2449 df_model : float 

2450 See model definition. 

2451 fitted_values : ndarray 

2452 Linear predictor XB. 

2453 llf : float 

2454 Value of the loglikelihood 

2455 llnull : float 

2456 Value of the constant-only loglikelihood 

2457 llr : float 

2458 Likelihood ratio chi-squared statistic; -2*(`llnull` - `llf`) 

2459 llr_pvalue : float 

2460 The chi-squared probability of getting a log-likelihood ratio 

2461 statistic greater than llr. llr has a chi-squared distribution 

2462 with degrees of freedom `df_model`. 

2463 prsquared : float 

2464 McFadden's pseudo-R-squared. 1 - (`llf`/`llnull`) 

2465 """ 

2466 

2467 def __init__(self, model, mlefit): 

2468 self.model = model 

2469 self.endog = model.endog 

2470 self.exog = model.exog 

2471 self.nobs = model.endog.shape[0] 

2472 

2473 # TODO: possibly move to model.fit() 

2474 # and outsource together with patching names 

2475 if hasattr(model, 'df_model'): 

2476 self.df_model = model.df_model 

2477 else: 

2478 self.df_model = len(mlefit.params) 

2479 # retrofitting the model, used in t_test TODO: check design 

2480 self.model.df_model = self.df_model 

2481 

2482 if hasattr(model, 'df_resid'): 

2483 self.df_resid = model.df_resid 

2484 else: 

2485 self.df_resid = self.endog.shape[0] - self.df_model 

2486 # retrofitting the model, used in t_test TODO: check design 

2487 self.model.df_resid = self.df_resid 

2488 

2489 self._cache = {} 

2490 self.__dict__.update(mlefit.__dict__) 

2491 

2492 def summary(self, yname=None, xname=None, title=None, alpha=.05): 

2493 """Summarize the Regression Results 

2494 

2495 Parameters 

2496 ---------- 

2497 yname : str, optional 

2498 Default is `y` 

2499 xname : list[str], optional 

2500 Names for the exogenous variables, default is "var_xx". 

2501 Must match the number of parameters in the model 

2502 title : str, optional 

2503 Title for the top table. If not None, then this replaces the 

2504 default title 

2505 alpha : float 

2506 significance level for the confidence intervals 

2507 

2508 Returns 

2509 ------- 

2510 smry : Summary instance 

2511 this holds the summary tables and text, which can be printed or 

2512 converted to various output formats. 

2513 

2514 See Also 

2515 -------- 

2516 statsmodels.iolib.summary.Summary : class to hold summary results 

2517 """ 

2518 

2519 top_left = [('Dep. Variable:', None), 

2520 ('Model:', None), 

2521 ('Method:', ['Maximum Likelihood']), 

2522 ('Date:', None), 

2523 ('Time:', None), 

2524 ('No. Observations:', None), 

2525 ('Df Residuals:', None), 

2526 ('Df Model:', None), 

2527 ] 

2528 

2529 top_right = [('Log-Likelihood:', None), 

2530 ('AIC:', ["%#8.4g" % self.aic]), 

2531 ('BIC:', ["%#8.4g" % self.bic]) 

2532 ] 

2533 

2534 if title is None: 

2535 title = self.model.__class__.__name__ + ' ' + "Results" 

2536 

2537 # create summary table instance 

2538 from statsmodels.iolib.summary import Summary 

2539 smry = Summary() 

2540 smry.add_table_2cols(self, gleft=top_left, gright=top_right, 

2541 yname=yname, xname=xname, title=title) 

2542 smry.add_table_params(self, yname=yname, xname=xname, alpha=alpha, 

2543 use_t=self.use_t) 

2544 

2545 return smry