GitHub PullRequestReviewCommentEvent at 2014-07-31 06:50:13
Repository: statsmodels
Actor: josef-pkt
Payload: {"action":"created","comment":{"url":"https://api.github.com/repos/statsmodels/statsmodels/pulls/comments/15628184","id":15628184,"diff_hunk":"@@ -0,0 +1,777 @@\n+\"\"\"\n+This module implements the Multiple Imputation through Chained Equations (MICE)\n+approach to handling missing data. This approach has 3 general steps:\n+\n+1) Simulate observations using a user specified conditional model.\n+2) Fit the model of interest to a complete, simulated dataset.\n+3) Repeat N times and combine the N models according to Rubin's Rules.\n+\n+Imputer instances, for imputing a single missing variable,\n+are specified with a (statsmodels) conditional model\n+(default is OLS with all other variables). A MICE instance is specified with\n+a model of interest together with its corresponding formulae. The results are\n+combined using the `combine` method.\n+\n+Reference for Rubin's Rules and Multiple Imputation:\n+\n+J L Schafer: \"Multiple Imputation: A Primer\", Stat Methods Med Res, 1999.\n+\n+Reference for Gaussian Approximation to the Posterior:\n+\n+T E Raghunathan et al.: \"A Multivariate Technique for Multiply\n+Imputing Missing Values Using a Sequence of Regression Models\",\n+Survey Methodology, 2001.\n+\n+Reference for Predictive Mean Matching:\n+\n+SAS Institute: \"Predictive Mean Matching Method for Monotone Missing Data\",\n+SAS 9.2 User's Guide, 2014.\n+\n+Reference for MICE Design in R package mi:\n+\n+A Gelman et al.: \"Multiple Imputation with Diagnostics (mi) in R: Opening\n+Windows into the Black Box\", Journal of Statistical Software, 2009.\n+\n+\"\"\"\n+#TODO: Add reference http://biomet.oxfordjournals.org/content/86/4/948.full.pdf\n+import operator\n+import pandas as pd\n+import numpy as np\n+import patsy\n+import statsmodels.api as sm\n+import random\n+import statsmodels\n+from statsmodels.tools.decorators import cache_readonly\n+from scipy import stats\n+import copy\n+import sys\n+\n+class ImputedData(object):\n+    __doc__= \"\"\"\n+    Stores missing data information and supports functionality for inserting\n+    values in missing data slots.\n+\n+    Can create Imputers directly via new_imputer method. By default, imputers\n+    are created for each variable using OLS using all other variables as \n+    predictors.\n+\n+    %(params)s\n+    data : array-like object\n+        Needs to support transformation to pandas dataframe. Missing value\n+        encoding is handled by pandas DataFrame class.\n+    go : boolean\n+        Governs whether to create Imputer objects automatically. If False, user \n+        must specify their own Imputers manually.\n+    method : string\n+        May take on values \"pmm\", \"gaussian\", or \"bootstrap\". Determines \n+        imputation method if go is True.\n+\n+    **Attributes**\n+\n+    data : pandas dataframe\n+        Dataset with missing values. After recording missing data information,\n+        simple column-wise means are filled into the missing values.\n+    columns : dictionary\n+        Stores indices of missing data.\n+    implist : MICE Imputer object\n+        One Imputer per variable to be imputed. See mice.Imputer.\n+\n+    \"\"\"\n+    def __init__(self, data, go=True, method=\"pmm\"):\n+        # may not need to make copies\n+        self.data = pd.DataFrame(data)\n+        # Drop observations where all variables are missing.\n+        self.data = self.data.dropna(how='all')\n+        self.columns = {}\n+        self.implist = []        \n+        for c in self.data.columns:\n+            self.columns[c] = MissingDataInfo(self.data[c])\n+            if go:\n+                self.new_imputer(c, method=method)\n+        # Fill missing values with column-wise mean.\n+        self.data = self.data.fillna(self.data.mean())\n+\n+    def new_imputer(self, endog_name, method=\"gaussian\", k_pmm=1, formula=None, \n+                    model_class=None, init_args={}, fit_args={}, \n+                    return_class=None, scale_method=\"fix\", scale_value=None, \n+                    transform=None, inv_transform=None):\n+        \"\"\"\n+        Create Imputer instance from our ImputedData instance.\n+        \n+        Adds imputer to self.implist.\n+\n+        Parameters\n+        ----------\n+        endog_name : string\n+            Name of the variable to be imputed.\n+        method : string\n+            May take on values \"pmm\", \"gaussian\", or \"bootstrap\". Determines \n+            imputation method, see mice.Imputer.       \n+        k_pmm : int\n+            Determines number of observations from which to draw imputation\n+            when using predictive mean matching. See mice.Imputer.\n+        formula : string\n+            Conditional formula for imputation. Defaults to model with main\n+            effects for all other variables in dataset.\n+        model_class : statsmodels model\n+            Conditional model for imputation. Defaults to OLS.\n+        return_class : scipy.random instance\n+            Controls the scale/location family to use as the asymptotic\n+            distribution of the imputed variable. For use with method\n+            impute_asymptotic_bayes.\n+        scale_method : string\n+            Governs the type of perturbation given to the scale parameter.\n+        scale_value : float\n+            Fixed value of scale parameter to use in simulation of data.\n+        transform : Numpy instance\n+            Transformation to apply to endogenous variable during imputation.\n+            Variable will be transformed back via inv_transform for the \n+            analysis model.\n+        inv_transform : Numpy instance\n+            Applied after imputation to recover original functional \n+            form of endogenous variable.\n+\n+        See Also\n+        --------\n+        mice.Imputer\n+        \"\"\"\n+        if model_class is None:\n+            model_class = sm.OLS\n+        if formula is None:\n+            formula = endog_name + \" ~ \" + \" + \".join(\n+                            [x for x in self.data.columns if x != endog_name])\n+        self.implist.append(Imputer(formula, model_class, self, method=method, \n+                                    k_pmm=k_pmm, init_args=init_args, \n+                                    fit_args=fit_args, \n+                                    scale_method=scale_method, \n+                                    scale_value=scale_value, \n+                                    transform=transform, \n+                                    inv_transform=inv_transform))\n+\n+    def store_changes(self, col, vals):\n+        \"\"\"\n+        Fill in dataset with imputed values.\n+\n+        Parameters\n+        ----------\n+        col : string\n+            Name of variable to be filled in.\n+        vals : array\n+            Array of imputed values to use in filling in missing values.\n+        \"\"\"\n+        \n+        ix = self.columns[col].ix_miss\n+        if len(ix) > 0:\n+            self.data[col].iloc[ix] = vals\n+\n+    def get_data_from_formula(self, formula):\n+        \"\"\"\n+        Use formula to construct endog and exog split by missing status.\n+\n+        Called by Imputer before fitting model.\n+\n+        Parameters\n+        ----------\n+        formula : string\n+            Patsy formula for the Imputer object's conditional model.\n+\n+        Returns\n+        -------\n+        endog_obs : DataFrame\n+            Observations of the variable to be imputed.\n+        exog_obs : DataFrame\n+            Observations of the predictors where the variable to be Imputed is\n+            observed.\n+        exog_miss : DataFrame\n+            Observations of the predictors where the variable to be Imputed is\n+            missing.\n+        \"\"\"\n+        dmat = patsy.dmatrices(formula, self.data, return_type=\"dataframe\")\n+        exog = dmat[1]\n+        endog = dmat[0]\n+        if len(endog.design_info.term_names) > 1:\n+            endog_name = tuple(endog.design_info.term_names)\n+        else:\n+            endog_name = endog.design_info.term_names[0]\n+        endog_obs = endog.iloc[self.columns[endog_name].ix_obs]\n+        exog_obs = exog.iloc[self.columns[endog_name].ix_obs]\n+        exog_miss = exog.iloc[self.columns[endog_name].ix_miss]\n+        return endog_obs, exog_obs, exog_miss\n+\n+class Imputer(object):\n+\n+    __doc__= \"\"\"\n+    Initializes object that imputes values for a single variable\n+    using a given formula.\n+\n+    %(params)s\n+\n+    formula : string\n+        Conditional formula used for imputation.\n+    model_class : statsmodels model\n+        Conditional model used for imputation.\n+    data : ImputedData object\n+        See mice.ImputedData\n+    method : string\n+        May take on values \"pmm\", \"gaussian\", or \"bootstrap\". Determines \n+        imputation method, see methods impute_pmm, impute_asymptotic_bayes, and \n+        impute_bootstrap.       \n+    k_pmm : int\n+        Determines number of observations from which to draw imputation\n+        when using predictive mean matching. See method impute_pmm.\n+    return_class : scipy.random class\n+        Controls the scale/location family to use as the asymptotic\n+        distribution of the imputed variable. For use in method\n+        impute_asymptotic_bayes.\n+    scale_method : string\n+        Governs the type of perturbation given to the scale parameter. See \n+        method perturb_params.\n+    scale_value : float\n+        Fixed value of scale parameter to use in simulation of data. See \n+        method perturb_params.\n+    transform : Numpy instance\n+        Transformation to apply to endogenous variable during imputation.\n+        Variable will be transformed back via inv_transform for the \n+        analysis model.\n+    inv_transform : Numpy instance\n+        Applied after imputation to recover original functional form of \n+        endogenous variable.\n+    %(extra_params)\n+\n+    **Attributes**\n+\n+    endog_name : string\n+        Name of variable to be imputed.\n+    num_missing : int\n+        Number of missing values.\n+\n+    Note: all params are saved as attributes.\n+\n+    \"\"\"\n+    def __init__(self, formula, model_class, data,\n+            method=\"gaussian\", k_pmm=1, init_args={}, fit_args={},\n+                 return_class=None, scale_method=\"fix\", scale_value=None, \n+                 transform=None, inv_transform=None):\n+        self.data = data\n+        self.formula = formula\n+        self.model_class = model_class\n+        self.init_args = init_args\n+        self.fit_args = fit_args\n+        self.endog_name = str(self.formula.split(\"~\")[0].strip())\n+        self.num_missing = len(self.data.columns[self.endog_name].ix_miss)\n+        self.return_class = return_class\n+        self.scale_method = scale_method\n+        self.scale_value = scale_value\n+        self.method = method\n+        self.k_pmm = k_pmm\n+        self.transform = transform\n+        self.inv_transform = inv_transform\n+\n+    def perturb_params(self, mdf):\n+        \"\"\"\n+        Perturbs the model's scale and fit parameters.\n+\n+        The user may choose to fix the scale at the estimated quantity or\n+        perturb it using a chi square approximation.\n+\n+        Bootstrap perturbation upcoming.\n+\n+        Parameters\n+        ----------\n+        mdf : statsmodels fitted model\n+            Passed from Imputer object.\n+\n+        Returns\n+        -------\n+        params : array\n+            Perturbed model parameters.\n+        scale_per : float\n+            Perturbed nuisance parameter.\n+        \"\"\"\n+        # TODO: switch to scipy\n+        params = mdf.params.copy()\n+        covmat = mdf.cov_params()\n+        covmat_sqrt = np.linalg.cholesky(covmat)\n+        if self.scale_method == \"fix\":\n+            if self.scale_value is None:\n+                scale_per = 1.\n+            else:\n+                scale_per = self.scale_value\n+        elif self.scale_method == \"perturb_chi2\":\n+            u = np.random.chisquare(float(mdf.df_resid))\n+            scale_per = float(mdf.df_resid) / float(u)\n+        elif self.scale_method == \"perturb_boot\":\n+            raise NotImplementedError\n+        p = len(params)\n+        params += np.dot(covmat_sqrt,\n+                         np.random.normal(0, np.sqrt(mdf.scale * scale_per), \n+                                          p))\n+        return params, scale_per\n+\n+    def impute_asymptotic_bayes(self):\n+        \"\"\"\n+        Use Gaussian approximation to posterior distribution to simulate data.\n+\n+        Fills in missing values of input data. User may also choose a different\n+        approximating location/scale family by specifying return_class in\n+        Imputer initialization.\n+        \"\"\"\n+        endog_obs, exog_obs, exog_miss = self.data.get_data_from_formula(\n+                                                                self.formula)\n+        if self.transform is not None and self.inv_transform is not None:\n+            endog_obs = self.transform(endog_obs)\n+        md = self.model_class(endog_obs, exog_obs, **self.init_args)\n+        mdf = md.fit(**self.fit_args)\n+        params, scale_per = self.perturb_params(mdf)\n+        new_rv = md.get_distribution(params=params, exog=exog_miss,\n+                                     model_class=self.return_class,\n+                                     scale=np.sqrt(scale_per * mdf.scale))\n+        new_endog = new_rv.rvs(size=len(exog_miss))\n+        if self.transform is not None and self.inv_transform is not None:\n+            new_endog = self.inv_transform(new_endog)\n+        self.data.store_changes(self.endog_name, new_endog)\n+\n+    def impute_pmm(self, pmm_neighbors=1):\n+        \"\"\"\n+        Use predictive mean matching to simulate data.\n+\n+        Fills in missing values of input data. Predictive mean matching picks\n+        an observation randomly from the observed value of the k-nearest\n+        predictions of the endogenous variable. Naturally, the candidate\n+        neighbors must have observed endogenous values.\n+\n+        Parameters\n+        ----------\n+        pmm_neighbors : int\n+            Number of neighbors in prediction space to select imputations from.\n+            Defaults to 1 (select closest neighbor).\n+        \"\"\"\n+        endog_obs, exog_obs, exog_miss = self.data.get_data_from_formula(\n+                                                                self.formula)\n+        if self.transform is not None and self.inv_transform is not None:\n+            endog_obs = self.transform(endog_obs)                                                                \n+        md = self.model_class(endog_obs, exog_obs, **self.init_args)\n+        mdf = md.fit(**self.fit_args)\n+        params, scale_per = self.perturb_params(mdf)\n+        # Predict imputed variable for both missing and nonmissing observations\n+        pendog_obs = md.predict(params, exog_obs)\n+        pendog_miss = md.predict(params, exog_miss)\n+        \n+#        imputed_miss = []\n+#        for val in pendog_miss:\n+#            dist = abs(pendog_obs - val)\n+#            dist_ind = np.argsort(dist, axis=0)\n+#            imputed_miss.append(random.choice(np.array(\n+#            endog_obs)[dist_ind[0:pmm_neighbors][0]]))\n+        \n+        ii = np.argsort(pendog_obs, axis=0)\n+        pendog_obs = pendog_obs[ii]\n+        oendog = endog_obs.iloc[ii,:]\n+        # Get indices of predicted endogs of nonmissing observations that are\n+        # close to those of missing observations\n+        ix = np.searchsorted(pendog_obs, pendog_miss)\n+        np.clip(ix, 0, len(pendog_obs) - 1, out=ix)\n+        ix_list = []\n+        for i in range(len(pendog_miss)):\n+            k = 0\n+            count_low = 0\n+            count_high = 0\n+#            ix_list.append([])\n+            upper = pendog_obs[ix[i]]\n+            lower = pendog_obs[ix[i] - 1]\n+            target = pendog_miss[i]\n+            ixs = []\n+#            ix_list[len(ix_list) - 1]\n+            limit_low = False\n+            limit_high = False\n+            while k < pmm_neighbors:\n+                if limit_low:\n+                    count_high += 1\n+                    ixs.append(ix[i] - 1 + count_high)\n+                elif limit_high:\n+                    count_low += 1\n+                    ixs.append(ix[i] - count_low)\n+                elif abs(target - upper) >= abs(target - lower):\n+                    count_low += 1\n+                    if ix[i] - 1 - count_low >= 0:\n+                        ixs.append(ix[i] - count_low)\n+                        target = copy.copy(lower)\n+                        lower = pendog_obs[ix[i] - 1 - count_low]\n+                    else:\n+                        ixs.append(ix[i] - 1 + count_high)\n+                        limit_low = True\n+                elif abs(target - upper) < abs(target - lower):\n+                    count_high += 1\n+                    if ix[i] + count_high < len(pendog_obs):\n+                        ixs.append(ix[i] - 1 + count_high)\n+                        target = copy.copy(upper)\n+                        upper = pendog_obs[ix[i] + count_high]\n+                    else:\n+                        ixs.append(ix[i] - count_low)\n+                        limit_high = True\n+                k += 1                    \n+            ixs = np.clip(ixs, 0, len(oendog) - 1)\n+            ix_list.append(random.choice(ixs))\n+        # Select a random draw of observed endogenous variable from a set\n+        # window around the closest predicted value\n+#        ix += np.random.randint(int(-pmm_neighbors / 2.), int(pmm_neighbors / 2.) + 1, size=len(ix))\n+#        np.clip(ix, 0, len(oendog), out=ix)\n+#        imputed_miss = np.array(oendog.iloc[ix,:])\n+        ix_list = np.squeeze(ix_list)\n+        imputed_miss = np.array(oendog.iloc[ix_list,:])\n+        if self.transform is not None and self.inv_transform is not None:\n+            imputed_miss = self.inv_transform(imputed_miss)        \n+        self.data.store_changes(self.endog_name, imputed_miss)\n+\n+    def impute_bootstrap(self):\n+        raise NotImplementedError\n+\n+class ImputerChain(object):\n+    __doc__= \"\"\"\n+    An iterator that returns imputed data sets produced using the MICE\n+    (multiple imputation by chained equations) procedure.\n+\n+    This class does imputation and returns the imputed data sets, it does not\n+    fit the analysis model. See the \"next\" method for details.\n+\n+    %(params)s\n+\n+    imputer_list : list\n+        List of Imputer objects, one for each variable to be imputed.\n+\n+    **Attributes**\n+\n+    data : pandas DataFrame\n+        Underlying data to be modified. Root copy of data is stored in original\n+        ImputedData object.\n+\n+    Note: All imputers must refer to the same data object. See mice.MICE.run\n+    for iterator call.\n+\n+    \"\"\"\n+    def __init__(self, imputer_list):\n+        self.imputer_list = imputer_list\n+        # Impute variable with least missing observations first\n+        self.imputer_list.sort(key=operator.attrgetter('num_missing'))\n+        # All imputers must refer to the same data object\n+        self.data = imputer_list[0].data.data\n+\n+    def __iter__(self):\n+        return self\n+\n+    def next(self):\n+        \"\"\"\n+        Makes this class an iterator that returns imputed datasets after\n+        cycling through all contained imputers.\n+\n+        Returned datsets are not saved unless specified in the iterator call.\n+\n+        Returns\n+        -------\n+        data : pandas DataFrame\n+            Dataset with imputed values saved after invoking each Imputer\n+            object in imputer_list.\n+        \"\"\"\n+        for im in self.imputer_list:\n+            if im.method==\"gaussian\":\n+                im.impute_asymptotic_bayes()\n+            elif im.method == \"pmm\":\n+                im.impute_pmm(im.k_pmm)\n+            elif im.method == \"bootstrap\":\n+                im.impute_bootstrap()\n+        return self.data\n+\n+class AnalysisChain(object):\n+    __doc__= \"\"\"\n+    An iterator that returns the fitted model of interest given a MICE imputed\n+    dataset.\n+\n+    Datasets to be used for analysis are chosen after an initial burnin period\n+    where no imputed data is used and also after skipping a set number of\n+    imputations for each iteration. See the \"next\" method for details.\n+\n+    Note: See mice.MICE.run for iterator call.\n+\n+    \"\"\"\n+\n+    def __init__(self, imputer_chain, analysis_formula, analysis_class,\n+                 skipnum=10, save=False, init_args={}, fit_args={}):\n+        self.imputer_chain = imputer_chain\n+        self.analysis_formula = analysis_formula\n+        self.analysis_class = analysis_class\n+        self.init_args = init_args\n+        self.fit_args = fit_args\n+        self.skipnum = skipnum\n+#        self.burnin = burnin\n+        self.save = save\n+        self.iter = 0\n+#        for b in range(self.burnin):\n+#            self.imputer_chain.next()\n+\n+    def __iter__(self):\n+        return self\n+\n+    def next(self):\n+        \"\"\"\n+        Makes this class an iterator that returns the fitted analysis model.\n+\n+        Handles skipping of imputation iterations, burnin period of\n+        unconsidered imputation iterations, and whether or not to save the\n+        datasets to which an analysis model is fit.\n+\n+        Returns\n+        -------\n+        mdf : statsmodels fitted model\n+            Fitted model of interest on imputed dataset that has passed all\n+            skip and burnin criteria\n+\n+        Note: If save option is True, imputed datasets are saved in the format\n+        \"mice_'iteration number'.csv\"\n+        \"\"\"\n+        #TODO: add transform here instead of in Imputer\n+        for i in range(self.skipnum):\n+            data = self.imputer_chain.next()\n+            print i\n+        md = self.analysis_class.from_formula(self.analysis_formula,\n+                                              data, **self.init_args)\n+        mdf = md.fit(**self.fit_args)\n+        if self.save:\n+            fname = \"%s_%d.csv\" % ('mice_', self.iter)\n+            data.to_csv(fname, index=False)\n+            self.iter += 1\n+        return mdf\n+\n+class MICE(object):\n+    __doc__= \"\"\"\n+    Fits the analysis model to each imputed dataset and combines the\n+    results using Rubin's rule.\n+\n+    Calls mice.Imputer_Chain and mice.AnalysisChain\n+    to handle imputation and fitting of analysis models to the correct imputed\n+    datasets, respectively.\n+\n+    %(params)s\n+\n+    analysis_formula : string\n+        Formula for model of interest to be fitted.\n+    analysis_class : statsmodels model\n+        Statsmodels model of interest.\n+    impdata : mice.ImputedData instance\n+        ImputedData object which contains an implist that is fully populated\n+        by the desired Imputer objects.\n+    %(extra_params)s\n+\n+    **Attributes**\n+\n+    N : int\n+        Total number of observations.\n+    modlist : list\n+        Contains analysis model fitted on each imputed dataset.\n+    df : list\n+        Contains degees of freedom for each parameter in the analysis model.\n+    fmi : list\n+        Contains fraction of missing information for each parameter in the \n+        analysis model.\n+    exog_names : list\n+        Contains names of exogenous variables in the analysis model.\n+    endog_name : string\n+        Name of endogenous model in the analysis model.\n+\n+    Examples\n+    --------\n+    >>>impdata = mice.ImputedData(data, go=False)\n+    >>>impdata.new_imputer(\"x2\", method=\"pmm\", k_pmm=20)\n+    >>>impdata.new_imputer(\"x1\", method=\"pmm\", model_class=sm.Logit, k_pmm=20)\n+    >>>impcomb = mice.MICE(\"x1 ~ x2\", sm.Logit,impdata)\n+    >>>impcomb.run(20, 10)\n+    >>>p1 = impcomb.combine()\n+    >>>print p1.summary()\n+\n+                               Results: MICE\n+==========================================================================\n+Method:              MICE              Dependent variable:             x1 \n+Model:               Logit             Sample size:                    978\n+--------------------------------------------------------------------------\n+           Coef.  Std.Err.    t     P>|t|   [0.025  0.975]    Df     FMI  \n+--------------------------------------------------------------------------\n+Intercept  2.2639   0.2182  10.3754 0.0000  1.8322  2.6956 130.0876 0.3412\n+x2        -3.8863   0.3304 -11.7609 0.5000 -4.5393 -3.2332 145.9336 0.3187\n+==========================================================================\n+    \"\"\"\n+    def __init__(self, analysis_formula, analysis_class, impdata,\n+                 init_args={}, fit_args={}):\n+        self.imputer_list = impdata.implist\n+        self.analysis_formula = analysis_formula\n+        self.analysis_class = analysis_class\n+        self.init_args = init_args\n+        self.fit_args = fit_args\n+        self.N = len(impdata.implist[0].data.data)\n+\n+    def run(self, num_ds=20, skipnum=10, burnin=30, save=False):\n+        \"\"\"\n+        Generates analysis model results.\n+\n+        Parameters\n+        ----------\n+        num_ds : int\n+            Number of imputed datasets to fit.\n+        skipnum : int\n+            Number of imputed datasets to skip between imputed datasets that\n+            are used for analysis. This is done to give the conditional\n+            distribution time to settle down. The literature says that the\n+            MICE procedure converges in distribution fairly quickly;\n+            standard practice is to set this number to be around ten.\n+        burnin : int\n+            Number of iterations to throw away before ever starting the skipped\n+            datasets count.\n+        save : boolean\n+            Whether to save the imputed datasets chosen for analysis.\n+        \"\"\"\n+       \n+        self.num_ds = num_ds\n+        imp_chain = ImputerChain(self.imputer_list)\n+        for b in range(burnin):\n+            imp_chain.next()         \n+        analysis_chain = AnalysisChain(imp_chain, self.analysis_formula,\n+                                       self.analysis_class, skipnum,\n+                                       save, self.init_args, self.fit_args)\n+        self.mod_list = []\n+        for current_iter in range(num_ds):\n+            achain = copy.deepcopy(analysis_chain)\n+            model = achain.next()\n+            self.mod_list.append(model)\n+            if not hasattr(self, \"exog_names\"):\n+                self.exog_names = model.model.exog_names\n+                self.endog_names = model.model.endog_names\n+            print current_iter\n+\n+    def combine(self):\n+        \"\"\"\n+        Pools estimated parameters and covariance matrices of generated\n+        analysis models according to Rubin's Rule.\n+\n+        Returns\n+        -------\n+        rslt : MICEResults instance\n+            Contains pooled analysis model results with missing data summaries.\n+        \"\"\"\n+\n+        params_list = []\n+        cov_list = []\n+        scale_list = []\n+#        full_cov = []\n+        for md in self.mod_list:\n+            params_list.append(md.params)\n+            cov_list.append(np.asarray(md.normalized_cov_params))\n+            scale_list.append(md.scale)\n+        scale = np.mean(scale_list)\n+        params = np.mean(params_list, axis=0)\n+        # Get average of within-imputation covariances weighted by scale\n+        full_cov = np.asarray(cov_list) * np.asarray(scale_list)[:, np.newaxis, np.newaxis]\n+        within_g = np.mean(full_cov, axis=0)\n+        # Used MLE rather than method of moments between group covariance\n+        between_g = np.cov(np.array(params_list).T, bias=1)\n+        cov_params = within_g + (1 + 1. / float(self.num_ds)) * between_g\n+#        gamma = (1. + 1. / float(self.num_ds)) * np.trace(np.dot(between_g,np.linalg.inv(cov_params))) \n+        gamma = (1. + 1. / float(self.num_ds)) * np.divide(np.diag(between_g), np.diag(cov_params))                \n+        df_approx = (float(self.num_ds) - 1.) * np.square(np.divide(1. , gamma))\n+#        np.sum(np.square(1. + np.diag(within_g)/(np.diag(between_g)*(1+1/float(self.N)))))\n+        df_obs = (float(self.N) - float(len(params)) + 1.) / (float(self.N) - float(len(params)) + 3.) * (1. - gamma) * (float(self.N) - float(len(params)))\n+        self.df = np.divide(1. , (np.divide(1. , df_approx) + np.divide(1. , df_obs)))\n+        self.fmi = gamma\n+#        self.within = within_g\n+#        self.between = between_g\n+        rslt = MICEResults(self, params, cov_params / scale)\n+        rslt.scale = scale\n+        return rslt\n+\n+\n+class MICEResults(statsmodels.base.model.LikelihoodModelResults):\n+\n+    def __init__(self, model, params, normalized_cov_params):\n+\n+        super(MICEResults, self).__init__(model, params,\n+                                          normalized_cov_params)\n+\n+    def summary(self, title=None, alpha=.05):\n+        \"\"\"\n+        Summarize the results of running MICE (multiple imputation with chained equations).\n+\n+        Parameters\n+        -----------\n+        title : string, optional\n+            Title for the top table. If not None, then this replaces\n+            the default title\n+        alpha : float\n+            Significance level for the confidence intervals\n+\n+        Returns\n+        -------\n+        smry : Summary instance\n+            This holds the summary tables and text, which can be\n+            printed or converted to various output formats.\n+        \"\"\"\n+\n+        from statsmodels.iolib import summary2\n+        from statsmodels.compat.collections import OrderedDict\n+        smry = summary2.Summary()\n+        float_format = \"%8.3f\"\n+\n+        info = OrderedDict()\n+        info[\"Method:\"] = \"MICE\"\n+        info[\"Model:\"] = self.model.analysis_class.__name__\n+        info[\"Dependent variable:\"] = self.model.endog_names\n+        info[\"Sample size:\"] = \"%d\" % self.model.mod_list[0].model.exog.shape[0]\n+\n+\n+        smry.add_dict(info, align='l', float_format=float_format)\n+\n+        param = summary2.summary_params(self, alpha=alpha)\n+        param['P>|t|'] = stats.t.sf(np.abs(np.asarray(param['t'])), self.model.df) / 2.\n+        ci = np.asarray(stats.t.interval(1-alpha, self.model.df, loc=np.asarray(param['Coef.']), scale=np.asarray(param['Std.Err.'])))\n+        param['[' + str(alpha/2)] = ci[0]\n+        param[str(1-alpha/2) + ']'] = ci[1]    \n+        param['Df'] = self.model.df\n+        param['FMI'] = self.model.fmi\n+        numiss = [0]\n+        for value in self.model.exog_names:\n+            for x in self.model.imputer_list:\n+                if x.endog_name == value:\n+                    numiss.append(int(x.num_missing))\n+#\n+#            t = next((x for x in self.model.imputer_list if x.endog_name == value), None)\n+#            numiss.append(t.num_missing)\n+        \n+        param['#missing'] = numiss\n+        smry.add_df(param, float_format=float_format)\n+        smry.add_title(title=title, results=self)\n+\n+        return smry\n+class MissingDataInfo(object):\n+    __doc__=\"\"\"\n+    Contains all the missing data information from the passed-in data object.\n+\n+    A self.columns dictionary entry exists for each column/variable in the\n+    dataset.\n+\n+    %(params)s\n+\n+    data : pandas DataFrame with missing values.\n+\n+    **Attributes**\n+\n+    ix_miss : array\n+        Indices of missing values for a particular variable.\n+    ix_obs : array\n+        Indices of observed values for a particular variable.\n+\n+    \"\"\"\n+\n+    def __init__(self, data):\n+        null = pd.isnull(data)\n+        self.ix_obs = np.flatnonzero(~null)\n+        self.ix_miss = np.flatnonzero(null)\n+        if len(self.ix_obs) == 0:\n+            raise ValueError(\"Variable to be imputed has no observed values\")\n+            ","path":"statsmodels/sandbox/mice/mice.py","position":776,"original_position":776,"commit_id":"3f82d568104d650a1dfa3a383b6d3015d1337a0b","original_commit_id":"3f82d568104d650a1dfa3a383b6d3015d1337a0b","user":{"login":"josef-pkt","id":440735,"avatar_url":"https://avatars.githubusercontent.com/u/440735?v=1","gravatar_id":"ad13088a623822caf74e635a68a55eae","url":"https://api.github.com/users/josef-pkt","html_url":"https://github.com/josef-pkt","followers_url":"https://api.github.com/users/josef-pkt/followers","following_url":"https://api.github.com/users/josef-pkt/following{/other_user}","gists_url":"https://api.github.com/users/josef-pkt/gists{/gist_id}","starred_url":"https://api.github.com/users/josef-pkt/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/josef-pkt/subscriptions","organizations_url":"https://api.github.com/users/josef-pkt/orgs","repos_url":"https://api.github.com/users/josef-pkt/repos","events_url":"https://api.github.com/users/josef-pkt/events{/privacy}","received_events_url":"https://api.github.com/users/josef-pkt/received_events","type":"User","site_admin":false},"body":"Is there trailing whitespace in this file?\r\nCan be cleaned during rebase without requiring a messy cleanup commit.","created_at":"2014-07-31T06:50:13Z","updated_at":"2014-07-31T06:50:13Z","html_url":"https://github.com/statsmodels/statsmodels/pull/1757#discussion_r15628184","pull_request_url":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757","_links":{"self":{"href":"https://api.github.com/repos/statsmodels/statsmodels/pulls/comments/15628184"},"html":{"href":"https://github.com/statsmodels/statsmodels/pull/1757#discussion_r15628184"},"pull_request":{"href":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757"}}},"pull_request":{"url":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757","id":17037321,"html_url":"https://github.com/statsmodels/statsmodels/pull/1757","diff_url":"https://github.com/statsmodels/statsmodels/pull/1757.diff","patch_url":"https://github.com/statsmodels/statsmodels/pull/1757.patch","issue_url":"https://api.github.com/repos/statsmodels/statsmodels/issues/1757","number":1757,"state":"open","title":"Mice","user":{"login":"chengfulun","id":6836993,"avatar_url":"https://avatars.githubusercontent.com/u/6836993?v=1","gravatar_id":"4ab66fb3667f38e83ab4242762f77256","url":"https://api.github.com/users/chengfulun","html_url":"https://github.com/chengfulun","followers_url":"https://api.github.com/users/chengfulun/followers","following_url":"https://api.github.com/users/chengfulun/following{/other_user}","gists_url":"https://api.github.com/users/chengfulun/gists{/gist_id}","starred_url":"https://api.github.com/users/chengfulun/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/chengfulun/subscriptions","organizations_url":"https://api.github.com/users/chengfulun/orgs","repos_url":"https://api.github.com/users/chengfulun/repos","events_url":"https://api.github.com/users/chengfulun/events{/privacy}","received_events_url":"https://api.github.com/users/chengfulun/received_events","type":"User","site_admin":false},"body":"This is a WIP PR. I have implemented a statsmodels class that receives a dataset with missing values and imputes the missing values using chained equations. It can then fit an analysis model to the fully imputed data. The code gives sensible results for a real dataset but has not been rigorously tested against existing software packages. ","created_at":"2014-06-12T04:40:56Z","updated_at":"2014-07-31T06:50:13Z","closed_at":null,"merged_at":null,"merge_commit_sha":"9e2be097e1c8ecd987f5ed05cfe23bbab21f2810","assignee":null,"milestone":null,"commits_url":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757/commits","review_comments_url":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757/comments","review_comment_url":"https://api.github.com/repos/statsmodels/statsmodels/pulls/comments/{number}","comments_url":"https://api.github.com/repos/statsmodels/statsmodels/issues/1757/comments","statuses_url":"https://api.github.com/repos/statsmodels/statsmodels/statuses/3f82d568104d650a1dfa3a383b6d3015d1337a0b","head":{"label":"chengfulun:mice","ref":"mice","sha":"3f82d568104d650a1dfa3a383b6d3015d1337a0b","user":{"login":"chengfulun","id":6836993,"avatar_url":"https://avatars.githubusercontent.com/u/6836993?v=1","gravatar_id":"4ab66fb3667f38e83ab4242762f77256","url":"https://api.github.com/users/chengfulun","html_url":"https://github.com/chengfulun","followers_url":"https://api.github.com/users/chengfulun/followers","following_url":"https://api.github.com/users/chengfulun/following{/other_user}","gists_url":"https://api.github.com/users/chengfulun/gists{/gist_id}","starred_url":"https://api.github.com/users/chengfulun/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/chengfulun/subscriptions","organizations_url":"https://api.github.com/users/chengfulun/orgs","repos_url":"https://api.github.com/users/chengfulun/repos","events_url":"https://api.github.com/users/chengfulun/events{/privacy}","received_events_url":"https://api.github.com/users/chengfulun/received_events","type":"User","site_admin":false},"repo":{"id":17358700,"name":"statsmodels","full_name":"chengfulun/statsmodels","owner":{"login":"chengfulun","id":6836993,"avatar_url":"https://avatars.githubusercontent.com/u/6836993?v=1","gravatar_id":"4ab66fb3667f38e83ab4242762f77256","url":"https://api.github.com/users/chengfulun","html_url":"https://github.com/chengfulun","followers_url":"https://api.github.com/users/chengfulun/followers","following_url":"https://api.github.com/users/chengfulun/following{/other_user}","gists_url":"https://api.github.com/users/chengfulun/gists{/gist_id}","starred_url":"https://api.github.com/users/chengfulun/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/chengfulun/subscriptions","organizations_url":"https://api.github.com/users/chengfulun/orgs","repos_url":"https://api.github.com/users/chengfulun/repos","events_url":"https://api.github.com/users/chengfulun/events{/privacy}","received_events_url":"https://api.github.com/users/chengfulun/received_events","type":"User","site_admin":false},"private":false,"html_url":"https://github.com/chengfulun/statsmodels","description":"Statsmodels: statistical modeling and econometrics in Python","fork":true,"url":"https://api.github.com/repos/chengfulun/statsmodels","forks_url":"https://api.github.com/repos/chengfulun/statsmodels/forks","keys_url":"https://api.github.com/repos/chengfulun/statsmodels/keys{/key_id}","collaborators_url":"https://api.github.com/repos/chengfulun/statsmodels/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/chengfulun/statsmodels/teams","hooks_url":"https://api.github.com/repos/chengfulun/statsmodels/hooks","issue_events_url":"https://api.github.com/repos/chengfulun/statsmodels/issues/events{/number}","events_url":"https://api.github.com/repos/chengfulun/statsmodels/events","assignees_url":"https://api.github.com/repos/chengfulun/statsmodels/assignees{/user}","branches_url":"https://api.github.com/repos/chengfulun/statsmodels/branches{/branch}","tags_url":"https://api.github.com/repos/chengfulun/statsmodels/tags","blobs_url":"https://api.github.com/repos/chengfulun/statsmodels/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/chengfulun/statsmodels/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/chengfulun/statsmodels/git/refs{/sha}","trees_url":"https://api.github.com/repos/chengfulun/statsmodels/git/trees{/sha}","statuses_url":"https://api.github.com/repos/chengfulun/statsmodels/statuses/{sha}","languages_url":"https://api.github.com/repos/chengfulun/statsmodels/languages","stargazers_url":"https://api.github.com/repos/chengfulun/statsmodels/stargazers","contributors_url":"https://api.github.com/repos/chengfulun/statsmodels/contributors","subscribers_url":"https://api.github.com/repos/chengfulun/statsmodels/subscribers","subscription_url":"https://api.github.com/repos/chengfulun/statsmodels/subscription","commits_url":"https://api.github.com/repos/chengfulun/statsmodels/commits{/sha}","git_commits_url":"https://api.github.com/repos/chengfulun/statsmodels/git/commits{/sha}","comments_url":"https://api.github.com/repos/chengfulun/statsmodels/comments{/number}","issue_comment_url":"https://api.github.com/repos/chengfulun/statsmodels/issues/comments/{number}","contents_url":"https://api.github.com/repos/chengfulun/statsmodels/contents/{+path}","compare_url":"https://api.github.com/repos/chengfulun/statsmodels/compare/{base}...{head}","merges_url":"https://api.github.com/repos/chengfulun/statsmodels/merges","archive_url":"https://api.github.com/repos/chengfulun/statsmodels/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/chengfulun/statsmodels/downloads","issues_url":"https://api.github.com/repos/chengfulun/statsmodels/issues{/number}","pulls_url":"https://api.github.com/repos/chengfulun/statsmodels/pulls{/number}","milestones_url":"https://api.github.com/repos/chengfulun/statsmodels/milestones{/number}","notifications_url":"https://api.github.com/repos/chengfulun/statsmodels/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/chengfulun/statsmodels/labels{/name}","releases_url":"https://api.github.com/repos/chengfulun/statsmodels/releases{/id}","created_at":"2014-03-03T07:42:55Z","updated_at":"2014-06-24T04:33:18Z","pushed_at":"2014-07-31T05:07:20Z","git_url":"git://github.com/chengfulun/statsmodels.git","ssh_url":"git@github.com:chengfulun/statsmodels.git","clone_url":"https://github.com/chengfulun/statsmodels.git","svn_url":"https://github.com/chengfulun/statsmodels","homepage":"http://statsmodels.sourceforge.net/devel/","size":19776,"stargazers_count":0,"watchers_count":0,"language":"Python","has_issues":false,"has_downloads":true,"has_wiki":true,"forks_count":0,"mirror_url":null,"open_issues_count":0,"forks":0,"open_issues":0,"watchers":0,"default_branch":"master"}},"base":{"label":"statsmodels:master","ref":"master","sha":"9da5dc5366bd6e148fbc9f84e77a30075ed3da71","user":{"login":"statsmodels","id":717666,"avatar_url":"https://avatars.githubusercontent.com/u/717666?v=1","gravatar_id":"cb5606256d549aade845966b3aae30d9","url":"https://api.github.com/users/statsmodels","html_url":"https://github.com/statsmodels","followers_url":"https://api.github.com/users/statsmodels/followers","following_url":"https://api.github.com/users/statsmodels/following{/other_user}","gists_url":"https://api.github.com/users/statsmodels/gists{/gist_id}","starred_url":"https://api.github.com/users/statsmodels/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/statsmodels/subscriptions","organizations_url":"https://api.github.com/users/statsmodels/orgs","repos_url":"https://api.github.com/users/statsmodels/repos","events_url":"https://api.github.com/users/statsmodels/events{/privacy}","received_events_url":"https://api.github.com/users/statsmodels/received_events","type":"Organization","site_admin":false},"repo":{"id":1885237,"name":"statsmodels","full_name":"statsmodels/statsmodels","owner":{"login":"statsmodels","id":717666,"avatar_url":"https://avatars.githubusercontent.com/u/717666?v=1","gravatar_id":"cb5606256d549aade845966b3aae30d9","url":"https://api.github.com/users/statsmodels","html_url":"https://github.com/statsmodels","followers_url":"https://api.github.com/users/statsmodels/followers","following_url":"https://api.github.com/users/statsmodels/following{/other_user}","gists_url":"https://api.github.com/users/statsmodels/gists{/gist_id}","starred_url":"https://api.github.com/users/statsmodels/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/statsmodels/subscriptions","organizations_url":"https://api.github.com/users/statsmodels/orgs","repos_url":"https://api.github.com/users/statsmodels/repos","events_url":"https://api.github.com/users/statsmodels/events{/privacy}","received_events_url":"https://api.github.com/users/statsmodels/received_events","type":"Organization","site_admin":false},"private":false,"html_url":"https://github.com/statsmodels/statsmodels","description":"Statsmodels: statistical modeling and econometrics in Python","fork":false,"url":"https://api.github.com/repos/statsmodels/statsmodels","forks_url":"https://api.github.com/repos/statsmodels/statsmodels/forks","keys_url":"https://api.github.com/repos/statsmodels/statsmodels/keys{/key_id}","collaborators_url":"https://api.github.com/repos/statsmodels/statsmodels/collaborators{/collaborator}","teams_url":"https://api.github.com/repos/statsmodels/statsmodels/teams","hooks_url":"https://api.github.com/repos/statsmodels/statsmodels/hooks","issue_events_url":"https://api.github.com/repos/statsmodels/statsmodels/issues/events{/number}","events_url":"https://api.github.com/repos/statsmodels/statsmodels/events","assignees_url":"https://api.github.com/repos/statsmodels/statsmodels/assignees{/user}","branches_url":"https://api.github.com/repos/statsmodels/statsmodels/branches{/branch}","tags_url":"https://api.github.com/repos/statsmodels/statsmodels/tags","blobs_url":"https://api.github.com/repos/statsmodels/statsmodels/git/blobs{/sha}","git_tags_url":"https://api.github.com/repos/statsmodels/statsmodels/git/tags{/sha}","git_refs_url":"https://api.github.com/repos/statsmodels/statsmodels/git/refs{/sha}","trees_url":"https://api.github.com/repos/statsmodels/statsmodels/git/trees{/sha}","statuses_url":"https://api.github.com/repos/statsmodels/statsmodels/statuses/{sha}","languages_url":"https://api.github.com/repos/statsmodels/statsmodels/languages","stargazers_url":"https://api.github.com/repos/statsmodels/statsmodels/stargazers","contributors_url":"https://api.github.com/repos/statsmodels/statsmodels/contributors","subscribers_url":"https://api.github.com/repos/statsmodels/statsmodels/subscribers","subscription_url":"https://api.github.com/repos/statsmodels/statsmodels/subscription","commits_url":"https://api.github.com/repos/statsmodels/statsmodels/commits{/sha}","git_commits_url":"https://api.github.com/repos/statsmodels/statsmodels/git/commits{/sha}","comments_url":"https://api.github.com/repos/statsmodels/statsmodels/comments{/number}","issue_comment_url":"https://api.github.com/repos/statsmodels/statsmodels/issues/comments/{number}","contents_url":"https://api.github.com/repos/statsmodels/statsmodels/contents/{+path}","compare_url":"https://api.github.com/repos/statsmodels/statsmodels/compare/{base}...{head}","merges_url":"https://api.github.com/repos/statsmodels/statsmodels/merges","archive_url":"https://api.github.com/repos/statsmodels/statsmodels/{archive_format}{/ref}","downloads_url":"https://api.github.com/repos/statsmodels/statsmodels/downloads","issues_url":"https://api.github.com/repos/statsmodels/statsmodels/issues{/number}","pulls_url":"https://api.github.com/repos/statsmodels/statsmodels/pulls{/number}","milestones_url":"https://api.github.com/repos/statsmodels/statsmodels/milestones{/number}","notifications_url":"https://api.github.com/repos/statsmodels/statsmodels/notifications{?since,all,participating}","labels_url":"https://api.github.com/repos/statsmodels/statsmodels/labels{/name}","releases_url":"https://api.github.com/repos/statsmodels/statsmodels/releases{/id}","created_at":"2011-06-12T17:04:50Z","updated_at":"2014-07-30T22:44:05Z","pushed_at":"2014-07-29T15:43:41Z","git_url":"git://github.com/statsmodels/statsmodels.git","ssh_url":"git@github.com:statsmodels/statsmodels.git","clone_url":"https://github.com/statsmodels/statsmodels.git","svn_url":"https://github.com/statsmodels/statsmodels","homepage":"http://statsmodels.sourceforge.net/devel/","size":39807,"stargazers_count":705,"watchers_count":705,"language":"Python","has_issues":true,"has_downloads":true,"has_wiki":true,"forks_count":275,"mirror_url":null,"open_issues_count":753,"forks":275,"open_issues":753,"watchers":705,"default_branch":"master"}},"_links":{"self":{"href":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757"},"html":{"href":"https://github.com/statsmodels/statsmodels/pull/1757"},"issue":{"href":"https://api.github.com/repos/statsmodels/statsmodels/issues/1757"},"comments":{"href":"https://api.github.com/repos/statsmodels/statsmodels/issues/1757/comments"},"review_comments":{"href":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757/comments"},"review_comment":{"href":"https://api.github.com/repos/statsmodels/statsmodels/pulls/comments/{number}"},"commits":{"href":"https://api.github.com/repos/statsmodels/statsmodels/pulls/1757/commits"},"statuses":{"href":"https://api.github.com/repos/statsmodels/statsmodels/statuses/3f82d568104d650a1dfa3a383b6d3015d1337a0b"}}}}