{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "How to generate missing values in Python? \n", "\n", "**Aude Sportisse with the help of Marine Le Morvan and Boris Muzellec**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Missing values occur in many domains and most datasets contain missing values (due to non-responses, lost records, machine failures, dataset fusions, etc.). These missing values have to be considered before or during analyses of these datasets.\n", "\n", "Now, if you have a method that deals with missing values, for instance imputation or estimation with missing values, how can you assess the performance of your method on a given dataset? If the data already contains missing values, than this does not help you since you generally do not have a ground truth for these missing values. So you will have to simulate missing values, i.e. you remove values – which you therefore know to be the ground truth – to generate missing values.\n", "\n", "The mechanisms generating missing values can be various but usually they are classified into three main categories defined by (Rubin 1976): missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR). The first two are also qualified as ignorable missing values mechanisms, for instance in likelihood-based approaches to handle missing values, whereas the MNAR mechanism generates nonignorable missing values. In the following we will briefly introduce each mechanism (with the definitions used widely in the literature) and propose ways of simulations missing values under these three mechanism assumptions. For more precise definitions we refer to references in the bibliography on the [R-miss-tastic](https://rmisstastic.netlify.app/bibliography/) website." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Notations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's denote by $\\mathbf{X}\\in\\mathcal{X_1}\\times\\dots\\times\\mathcal{X_p}$ the complete observations. We assume that $\\mathbf{X}$ is a concatenation of $p$ columns $X_j\\in\\mathcal{X_j}$, $j\\in\\{1,\\dots,p\\}$, where $dim(\\mathcal{X_j})=n$ for all $j$. \n", "\n", "The data can be composed of quantitative and/or qualitative values, hence $\\mathcal{X_j}$ can be $\\mathbb{R}^n$, $\\mathbb{Z}^n$ or more generally $\\mathcal{S}^n$ for any discrete set $S$.\n", "\n", "Missing values are indicated as `NA` (not available) and we define an indicator matrix $\\mathbf{R}\\in\\{0,1\\}^{n\\times p}$ such that $R_{ij}=1$ if $X_{ij}$ is observed and $R_{ij}=0$ otherwise. We call this matrix $\\mathbf{R}$ the response (or missingness) pattern of the observations $\\mathbf{X}$. According to this pattern, we can partition the observations $\\mathbf{X}$ into observed and missing: $\\mathbf{X} = (\\mathbf{X}^{obs}, \\mathbf{X}^{mis})$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Definition of the mechanisms \n", "\n", "In order to define the different missing values mechanisms, both $\\mathbf{X}$ and $\\mathbf{R}$ are modeled as random variables with probability distributions $\\mathbb{P}_X$ and $\\mathbb{P}_R$ respectively. We parametrize the missingness distribution $\\mathbb{P}_R$ by a parameter $\\phi$.\n", "\n", "### MCAR \n", "\n", "The observations are said to be Missing Completely At Random (MCAR) if the probability that an observation is missing is independent of the variables and observations: the probability that an observation is missing does not depend on $(\\mathbf{X}^{obs},\\mathbf{X}^{mis})$. Formally this is:\n", "$$\\mathbb{P}_R(R\\,|\\, X^{obs}, X^{mis}; \\phi) = \\mathbb{P}_R(R) \\qquad \\forall \\, \\phi.$$\n", "\n", "### MAR\n", "\n", "The observations are said to be Missing At Random (MAR) if the probability that an observation is missing only depends on the observed data $\\mathbf{X}^{obs}$. Formally,\n", "\n", "$$\\mathbb{P}_R(R\\,|\\,X^{obs},X^{mis};\\phi)=\\mathbb{P}_R(R\\,|\\,X^{obs};\\phi) \\qquad \\forall \\,\\phi,\\, \\forall \\, X^{mis}.$$\n", "\n", "### MNAR\n", "\n", "The observations are said to be Missing Not At Random (MNAR) in all other cases, i.e. the missingness depends on the missing values and potentially also on the observed values.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Use of `produce_NA` with default settings\n", "\n", "For now, with the main function `produce_NA`, it is possible to generate missing values only for quantitative data which are complete.\n", "\n", "Missing values can be generated following one or more of the three main missing values mechanisms (see below for details).\n", "\n", "The function is widely based on the code of Boris Muzellec available [here](https://github.com/BorisMuzellec/MissingDataOT)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We generate a small example of observations $\\mathbf{X}$:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "!pip install wget \n", "\n", "import wget\n", "wget.download('https://raw.githubusercontent.com/BorisMuzellec/MissingDataOT/master/utils.py')\n", "\n", "import numpy as np\n", "import pandas as pd\n", "from utils import *\n", "import torch\n", "import seaborn as sns" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Fix the seed ------------------------------------------------------\n", "np.random.seed(0)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0123
02.3801724.4053631.0F
13.3791802.4008690.0D
21.3415824.9441510.0E
31.4501412.9228430.0D
41.2660930.6957262.0D
\n", "
" ], "text/plain": [ " 0 1 2 3\n", "0 2.380172 4.405363 1.0 F\n", "1 3.379180 2.400869 0.0 D\n", "2 1.341582 4.944151 0.0 E\n", "3 1.450141 2.922843 0.0 D\n", "4 1.266093 0.695726 2.0 D" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Sample data generation ------------------------------------------------------\n", "# Generate complete data\n", "\n", "n = 100\n", "p = 2\n", "\n", "mu_X = np.repeat(1, p)\n", "Sigma_X = np.array([[1, 1], [1, 4]])\n", "X_complete_cont = np.random.multivariate_normal(mu_X, Sigma_X, size=n)\n", "\n", "from scipy.stats import poisson\n", "\n", "lamb = 0.5\n", "X_complete_discr = poisson.rvs(lamb, size=n)\n", "X_complete_discr = np.expand_dims(X_complete_discr, axis=1)\n", "\n", "\n", "n_cat = 5\n", "X_complete_cat = np.random.binomial(n=n_cat, p=0.5, size=n)\n", "X_complete_cat = np.expand_dims(X_complete_cat, axis=1)\n", "\n", "X_complete = np.concatenate((X_complete_cont, X_complete_discr, X_complete_cat),axis=1)\n", "\n", "X_complete = pd.DataFrame(X_complete)\n", "X_complete.iloc[:, 3] = X_complete.iloc[:, 3].astype('category')\n", "X_complete.iloc[:, 3].cat.categories = [\"F\", \"E\", \"D\", \"C\", \"B\", \"A\"]\n", "X_complete.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Minimal set of arguments\n", "\n", "In order to generate missing values for given data, `produce_NA` requires the following arguments:\n", "\n", "* `X`: the initial data (can be only complete for now) as a matrix or data.frame\n", "\n", "* `p_miss`: proportion of missing values to generate for variables which will have missing values\n", "\n", "* `mecha`: one of \"MCAR\", \"MAR\", \"MNAR\" (default: “MCAR”)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Function produce_NA for generating missing values ------------------------------------------------------\n", "\n", "def produce_NA(X, p_miss, mecha=\"MCAR\", opt=None, p_obs=None, q=None):\n", " \"\"\"\n", " Generate missing values for specifics missing-data mechanism and proportion of missing values. \n", " \n", " Parameters\n", " ----------\n", " X : torch.DoubleTensor or np.ndarray, shape (n, d)\n", " Data for which missing values will be simulated.\n", " If a numpy array is provided, it will be converted to a pytorch tensor.\n", " p_miss : float\n", " Proportion of missing values to generate for variables which will have missing values.\n", " mecha : str, \n", " Indicates the missing-data mechanism to be used. \"MCAR\" by default, \"MAR\", \"MNAR\" or \"MNARsmask\"\n", " opt: str, \n", " For mecha = \"MNAR\", it indicates how the missing-data mechanism is generated: using a logistic regression (\"logistic\"), quantile censorship (\"quantile\") or logistic regression for generating a self-masked MNAR mechanism (\"selfmasked\").\n", " p_obs : float\n", " If mecha = \"MAR\", or mecha = \"MNAR\" with opt = \"logistic\" or \"quanti\", proportion of variables with *no* missing values that will be used for the logistic masking model.\n", " q : float\n", " If mecha = \"MNAR\" and opt = \"quanti\", quantile level at which the cuts should occur.\n", " \n", " Returns\n", " ----------\n", " A dictionnary containing:\n", " 'X_init': the initial data matrix.\n", " 'X_incomp': the data with the generated missing values.\n", " 'mask': a matrix indexing the generated missing values.s\n", " \"\"\"\n", " \n", " to_torch = torch.is_tensor(X) ## output a pytorch tensor, or a numpy array\n", " if not to_torch:\n", " X = X.astype(np.float32)\n", " X = torch.from_numpy(X)\n", " \n", " if mecha == \"MAR\":\n", " mask = MAR_mask(X, p_miss, p_obs).double()\n", " elif mecha == \"MNAR\" and opt == \"logistic\":\n", " mask = MNAR_mask_logistic(X, p_miss, p_obs).double()\n", " elif mecha == \"MNAR\" and opt == \"quantile\":\n", " mask = MNAR_mask_quantiles(X, p_miss, q, 1-p_obs).double()\n", " elif mecha == \"MNAR\" and opt == \"selfmasked\":\n", " mask = MNAR_self_mask_logistic(X, p_miss).double()\n", " else:\n", " mask = (torch.rand(X.shape) < p_miss).double()\n", " \n", " X_nas = X.clone()\n", " X_nas[mask.bool()] = np.nan\n", " \n", " return {'X_init': X.double(), 'X_incomp': X_nas.double(), 'mask': mask}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Value\n", "\n", "`produce_NA` returns a list containing three elements: \n", "\n", "- `X_init`: the initial data\n", "- `X_incomp`: the data with the newly generated missing values\n", "- `mask`: a matrix indexing the generated missing values\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percentage of newly generated missing values: 39.0 %\n" ] } ], "source": [ "# Minimal example for generating missing data ------------------------\n", "X_miss_mcar = produce_NA(X_complete_cont, p_miss=0.4, mecha=\"MCAR\")\n", "\n", "X_mcar = X_miss_mcar['X_incomp']\n", "R_mcar = X_miss_mcar['mask']\n", "\n", "print(\"Percentage of newly generated missing values: \", (R_mcar.sum()).numpy()/np.prod(R_mcar.size())*100, \" %\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the data as a color-encoded matrix, cells with missing values are masked (in white). " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAD7CAYAAADNT5fNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAeRklEQVR4nO3de5RU5Znv8e8jxqioQY3cETV6vEczchgNRlDwAsFrjBGDt4xDciI55mpQ1xxNspwhMRn1LJ0Z0BjveEEx3hVQQD3RIzhRo8Y5inJroJXoeI/SPOeP2ozVu6u79q79Vu1d1b/PWr2km6p3v7Kqnn7q3e/zPubuiIhIY22S9wRERHojBV8RkRwo+IqI5EDBV0QkBwq+IiI5UPAVEclBpuBrZkeZ2ctm9oqZTQs1KRGRVme17vM1sz7AfwCHAyuBp4FJ7v5iuOmJiLSmTTM8dyTwirsvBTCzW4BjgW6Dr5mpokNEEnF3yzrGJ28uTRRzPvP5XTJfK60swXcIsKLs+5XA38YfZGZTgCkAl31lL87cc2iGS0qr+dzMh/lg4e/ynoYUzJajz8x7CnWXJfhW+k3R5beMu88EZkIp8/3+Y1qVkM56wxtNcrKhI+8ZdCtL8F0JDCv7fijQ1tMT3v+XqRkuJ62o73ev4J1vH5n3NKRgtpnxUJiBOtaHGacOstxw25TSDbexwCpKN9xOcfcXeniO1nxFJJEQa74ft72QKOZsNnjv5lnzdff1ZjYVeAjoA1zTU+AFWHbAYbVeTlrU8CWPoJP1JM4sUCzcsCHMOHVQc+Zb08WU+YpIQkEy3xXPJst8h+3XPJlvLZThSJyZ8d7PTsl7GlIwW114c5iBWvSGm4hIsXmLLjuY2evAu0AHsN7dR1R5vFJfEUkkxLLDX199MlHM+ewXDmzKZYdD3f3NJA/UsoPEmRlLhh6T9zSkYA5YeXeYgQp8w62hyw7B7mBKSwn2RhOJK/CyQ9bg68DD0XLCjKiarZNO5cUH7cEZuw/JeElpJf2unc8Lu0zIexpSMHsvvT/MQAW+4ZZ1zXewu7eZWX9gLvA9d1/Uw+O17iAiiQRZ833p0WRrvnse2lxrvu7eFv233czmUDrprNvgu/oro7NcTlrQoMcW8s7Z4/OehhTMNlc+EGagFi0v7gts4u7vRn+eC/zc3R/s4TnKfEUkkSCZ73MPJct8v3hkU2W+A4A50U20TYGbewq8AG2jxmS4nLSiwU8s4I2jDsl7GlIwOzzY7QfoVNxbdM039cWU+YpIQiEy34/+eG+imLP5/hObKvNNrX2cMhzprP+8RXxw1y/znoYUzJbH/TTMQAXe56vMV0QKKUjmu+SuZJnvAccVL/M1s2uAiUC7u+8T/Ww74FZgJ+B14CR3f6vaWKpwkzgz480J+kQknX3+/jBrvnR8EmacOkjSOv5a4KjYz6YB8919N2B+9L2ISLFs2JDsKwdVM193X2RmO8V+fCwwJvrzdcACoOoijcqLpZJgWY5IXAuWFw9w99UA7r46qnCrqLy8+H9s/d85Ystda7yktKLj186iz2cG5z0NKZiOT3psB5lcgW+41X23Q7x78b+++3S9LylNJtgbTSSuBYPvWjMbFGW9g4D2JE/SASoSt/fS+3nvwpPznoYUzFY/uyXIOF7gG261Bt+7gdOB6dF/f5/kScFOKpKWEuqNJtJFM6/5mtksSjfXPm9mK4ELKQXd28zs74DlwNeTXExbzSTOzPS6kC56Q/fiJLsdJnXzV2PTXky7HaQSvS6kbpo58w3p+Z2/2sjLSRPY97X7lPlKF8p8A9v3tfsaeTlpEsp8pW6aOfPtprz4IuDvgTeih53v7lXvpinDkTit+UolwX4hrw97mLqZ9QEWA6vcfWKWsZJkvtcCVwDXx35+qbv/Os3FlOFIJXpdSN2Ez3zPAV4Ctsk6UK3lxTVZPuKwEMNIC9lx8SMs/eK4vKchBbPLc/PCDBRwzdfMhgJfBS4Gfph1vCQH63Rnqpk9Z2bXmNm23T3IzKaY2WIzW3zTG6syXE5EJCXfkOirPE5FX1MqjHYZcC4QJKInOs83ynzvLVvzHQC8Sal1/C+AQe7+rQTjaHFPRBIJcZ7vh3OmJ4o5Wxw/rcdrmdlEYIK7f9fMxgA/bsSabxfuvrZsUlcB9yZ8Xi2XkxZmZtww6Jt5T0MK5tTVN4UZKNya7yjgGDObAGwObGNmN7r75FoHrCn4bjzXIfr2eOBPCZ9Xy+WkxQV7o4nEBdrt4O7nAecBlGW+NQdeqL28eIyZ7U9p2eF14NtJLqaDdSRu76X36xORdBEsUSvwa0s93ESkkIKs+c66MNma76SfFa+HW0jvnnt8Iy8nTWDrX81R5itdqLw4sK1/NaeRl5MmoXsBUjdNXl48jFJ120BK+9tmuvvltXQwVoYjcWbG/P6JTiSVXmRs++1hBuroCDNOHSTJfNcDP3L3Z8xsa2CJmc0FzqDUwXi6mU2j1MG4xyaaynCkkmBvNJG4Zl52iLaUbWyW+a6ZvQQMoYYOxsp8Jc7M+Odh2ucrnf1wRaDthwUOvqnKi6NKty8BTxHrYAxU7GBcXrY3c+bMbLMVEUkjYXlxHhJvNTOzrYCFwMXufqeZve3u/cr+/i137/aMh+gxSn1FJJEQW80+mPmDRDFnyymXFnOrmZl9BrgDuMnd74x+nLqDsZYdJE7n+Uol2moGWOlf4bfAS+7+z2V/lbqDsW64SSV6XUjdNPluh1HAqcDzZvbH6GfnU0MH40f7n1jrPKVFHdo+m3v7n5z3NKRgJrbfEmagAme+Ki8WkUIKsuZ7+XeSrfme82/FXPMNRd2LJW7f1+5jynB9IpLOZi6bHWagAt9PUPdiyV2wN5pIXIGXHbKUF19Eyg7Gtw08JdtspeWctOZm7XaQLsLtdijuaytLeTGk7GB80pqba5mjtDjtdpC6aebdDj2UF6f27zseXcvTpIV9afk9/HJHlRdLZz9dHqa82Jt52aFcrLx4FKUOxqcBiyllx11ONYu6gFbqBCoChHujiXTR5MsOwH+VF98BfN/d3zGzf6XUuXhjB+PfAF06GLv7TGDjoQ7F/ZeQXJgZ8/qflPc0pGDGtd8WZqBmPs8XKpcX19rBWESkYZo58+2uvLiWDsa6sSKVBMtyROLWN/ENN7ovL56UtoOxthRJnJnx3E4T856GFMwXXw/0QbqZlx3c/XGgUsra457eSpT5SiXB3mgicc287BDSa/uNa+TlpAns/Ow8Xt7jyLynIQWz+58fCjJOU281M7PNgUXAZ6PHz3b3C2tpoLnzs/OyzldaUKg3mkgXTZ75/hU4zN3fi3Y9PG5mDwAnkLKBprrUStzY9tu5RWXnEnNyqGrYZg6+XrpL9l707WeiL6eGBprqUiuVBHujicQ1c3kxgJn1AZYAuwJXuvtTZtapgaaZVWygWU67HSTOzLhxkMqLpbPJq0OVFxc35iQKvu7eAexvZv2AOWa2T9ILqLxYqgn1RhPpotmD70bu/raZLQCOImEDTZUXS0/UQFMqUQNNwMx2AD6JAu8WwDjgl6iBpgSi14XUTZNnvoOA66J1302A29z9XjP7AykbaOpISYn70vJ7WHngoXlPQwpm6JOPhhmomYOvuz9H6RjJ+M/XAWPrMSkRkRC8o4mXHUL60vJ7Gnk5aRLBshyRuGbOfEPSjRWJMzPGDxuf9zSkYB5Y8UCQcZp6q1kP5cUXkbKBpm6sSCWh3mgiXTRz8KX78mJI2UDziYEn1DJHaWGj1tzJJcNUZCGd/WRFoL3fgZZ8u+vinmXMLOXFqY1ac2ctT5MWF+yNJhLj64PdcKvYxd3dX6x1wCzlxeNJ0ECznNZ8Jc7MeOfbOlJSOttmRqCT7gLF3h66uNccfC1NQNxYXgx8j9Ja75t82kBzkLt3aaAZKy8+oNaJikjv4u6ZbxK99fUxiQLctrcvSHytqIv7ImAfd3+ntpllKC8uX+vtqYFmeXnxM8OOVeornRyw8m59IpIuwpUXJ75e/AyamVHsij+uUxf3LFOruby4lgaaB6y8O8tcpUVpF4zUS9KtZrEzaCqq1MU9iyzlxTeogaZkZWasO+6QvKchBbP9XYvCDBRut0PFLu6ZxmxkQDQzRV8RSSTEmu+6r45OFHO2v29hj9cys4OBx4Dn+TSkV61t6Ikq3CRXOlJSKgm1FBWqc3wPXdxr1tDgKyLSUMU9Vyd58I3WfBcDq9x9Yi3di3VjRSrR60LqJVTmWw9pMt9zgJeAbaLvp5Gye7E+XkqcmfFo/xPznoYUzKHts4OM0/TB18yGAl8FLgZ+GP04dfdiZThSSag3mkicdxQ35iTNfC8DzgW2LvtZou7F5ZuXZ8yYwZQp6qUpnzIzzh5+Ut7TkIK5ctltQcYpcuZbdauZmU0EJrj7d81sDPDjaM33bXfvV/a4t9x92ypjad1BRBIJsdVs9cGHJoo5gx5/tOEpcpLMdxRwjJlNADYHtjGzG0nYvbjc3QNOzjZbaTnHrL1F9wKki6JtNauHtAfrjOHTzPcSYF3ZDbft3P3cKs/Xu0xEEgmR+a466LBEMWfIHx4pZObbnemk7F78/r9MzXA5aUV9v3uFMl/pQplv6Isp8xWRhEJkvstHjE0Uc3ZcPL+pMt/UVn9ldCMvJ01g0GMLlflKF+Ey3+bfaiYi0nRaIvhWKC++iJTdiwc9trDWeUoLU/GN1EuRP1RlKS+GlN2L9fFS4syMHwz/Rt7TkIK5dNmtQcZp+sy3m/Li1JThSCWh3mgicQHu2dVNlvJiSNC9uLy8+Nd77MZpQwdlmK60mv7zFnH7wFPynoYUzNfX3BxknI4Cn+2Qpbx4AAm6F8fG0rqDiCQSYqvZy3uMTxRzdv/zA4XcalaxvNjdJ298QE/di8tpzVfi1MlCKtFWM8DdzwPOg07lxZNr6V6sNV+pRK8LqZci/17Pss/3V+peLFmZGW2jxuQ9DSmYwU8sCDJOkTNflReLSCGFWPN9fuejE8WcfV+7p5BrvsEo85U4M+PJQcfnPQ0pmANXzwkyTpFDTkODr9b2pJJQbzSRuA3Nvs/XzF4H3gU6gPXuPqKW7sXLRxyWZa7SgnZc/Aiv7H1E3tOQgtn1hYeDjFPkIotNUjz2UHff391HRN9v7F68GzA/+l5EpDDck33lIdENtyjzHeHub5b97GVgTFkboQXuvnuVcQq8AiMiRRLihtviocclijkjVt5V2BtuDjwcBc8Z7j4TdS+WAMyMj/7jibynIQWz+X8bFWScjg1pPtw3VtLMd7C7t0UBdi7wPeBudS8WkXoJkfk+OfiERDHnwLY7i5n5untb9N92M5sDjKSG7sUv7TY+02Sl9ez5/x7gzOFfy3saUjC/W3ZHkHGKvNshycE6fYFN3P3d6M9zgZ8DY1H3YhGpkxCZ7xMDT0wUc0atmV3IzHcAMCfao7spcLO7P2hmT5Oye7GKLCTOzDhj+Al5T0MK5tpldwYZp8DNi1VeLCLFFCLzXTTw64liziFrbi9k5hvMA/3VLkY6G99+qz4RSRehqmHXF3jNt6HBd3y72sVIVyo7l3pxivvaylJefBEpuxe/f8XZtc9UWlLfqVfy4aLr856GFMwWh5wWZJwir/mmyXwPLa9wi6TqXtx36pUpLie9Rag3mkhc02e+obz3D1rzlc62+sWtfPj4jXlPQwpmi4MnV39QAiEzXzM7Crgc6ANc7e7Ts4yXtPZuY3nxkqhceKOpZvacmV1jZhWr28xsipktNrPF1zzzapa5ioik0oEl+qrGzPoAVwLjgb2ASWa2V5a5ZSkvfhl1LxaROgmx1eyegZMSxZyj18zq8VpmdhBwkbsfGX1/XjTHf6p1bjWXF7v7orKJJepe/JcTR9c6T2lR281eyAXDJ+U9DSmYi5fNCjLOhoRrvuUHgEVmRgeIbTQEWFH2/Urgb7PMrWrwrVBefATw81q6F283e2GWuUqLCvVGE4lL+lE7CrQze3hIpSie6ZN8lvLiG9S9WLIyM+b1PynvaUjBjGu/Lcg4AW+4rQSGlX0/FGjLMmDV4OvuS4H9Kvz81LQX02Z6qSTUG00kbkO4mPM0sJuZ7QysAk4GTskyoLoXS67MTK8L6SJUotYRZBRw9/VmNhV4iNJWs2vc/YUsY6p7seROrwuplw0BX1pRBW+PVbxpJC0v7gdcDexDaY33W5S2mqXqXvzWyWNqn6m0pG1vWcBfTtAuGOlsuzvD3JxPutshD0n3+V4HPObuV5vZZsCWwPnAX8oOU9/W3X9aZRx9vhSRRELs871x8OREMWdy243FO1LSzLYBDgHOAHD3j4GPzexYYEz0sOuABUCPwffR/ifWPlNpSYe2z9b+b+ki1LbUkMsOoSVZdtiF0sllvzOz/YAlwDnU0L1YpBLt/5Z6afZTzTYF/gb4nrs/ZWaXA9OSXiC2eVnLDtKJmfHafuPynoYUzM7PzgsyTkeBM98kB+usBFa6+1PR97MpBeO1UddiknYvFhFppA0Jv/KQpMhijZmtMLPd3f1lSl2LX4y+TgemR//9fbWxtKVIKgmV5YjENfuyA5ROMbsp2umwFDiTUtas7sWSiZlx3eBv5j0NKZjT224KMk6BW7glPtXsj8CICn81Ns3FlPlKJaHeaCJxrZD5BvH2aYc18nLSBPpd/wg3DFLmK52dujrML+RQ5cX10NDg2+/6Rxp5OWkSod5oInHNvs+3u/LiI0nZvVhrvhJnZrx9RqrVK+kF+l07P8g4rbDscDnwoLufWFZefCQpuxdrzVcqCfVGE4lr6uDbQ3lx6ou9uOv41M+R1rbXKw/oE5F0ESpRK/IrK0t5MZS6F58GLAZ+VOlUM5UXSzX6RCT1UuQ136qnmpnZCOBJYFRZefE7wBWk7F789JDji/yLSHIwsu0u7fOVLk5vuynIqWb/NDzZqWbnLWv8qWY1lxe7+1p373D3DcBVwMh6TVJEpBYb8ERfeai5vLiW7sUj2+7KNltpSSqykHpp6htukUrlxf9b3YslK/Vwk0p0wy3STXmxuhdLEHpdSL20QuYbhDIciTMz3p6ssnPprN+NYaph1xe4c1mSfb67U2qUudEuwP8CridlA01lOFJJqDeaSFxxQ2+yG24vA/sDmFkfYBUwh1I3i/llDTSnUaWHmzJfidOar1QSKlFrpWWHscCr7r6slgaaynylEr0upF7y2kaWRNrgezIwK/pzogaa5R4f8LWUl5NWd/DaO5T5Shfa7VAm2mZ2DHBemguovFiqUeYr9dIqyw7jgWfcfW30/dqNhRY9NdBU92LpidZ8pZJQv5A7Chxy0gTfSXy65ABwN2qgKQHodSH10vSZr5ltCRxO5yq26aRsoPn+b86qZY7Swvr+6Go+fuPVvKchBbPZDl8IMo43e+br7h8A28d+to6UDTRFRBqp6TPfUPr+6OpGXk6aRKgsRySulbaaZbLsAJWRSmfDlzzChTuekvc0pGB+tvzmIOMUN/RmKy/uR8oGmsOXqIxUugr1RhOJW1/g8JulvPhMUjbQ1JYiidNWM6kkXJFFcV9bWcqLU19MW4qkEr0upF5a6YZbeXkxJGigWe4/zzo85eWk1X3u6rkqO5cuDl57R5Bxipz5Vm2g+V8PLJUXtwF7u/taMxtAggaasfLiA4LMWkRaXogGmqfv9LVEAe661+/IdC0zuwQ4GvgYeBU4093f7uk5NZcXl5UZY2ZXAfdWepLKi6UnZsYNg9S9WDo7dXWYvn4djbufMBc4z93Xm9kvKZ2B0+MpjzWXF9fSQFNre1JJqDeaSFyj9vm6+8Nl3z4JnFjtOVnKi3+lBpqSlZnx4q7j856GFMxerzwQZJyka74VTl+cGX1qr8W36Lw9t/I1GxkQzQrcUElECiXEmu83hh+XKObcuuyuqtcys3nAwAp/dYG7/z56zAWUmg2f4FWCa0Mr3J4dPrGRl5MmsN+ye/ngqh/kPQ0pmC3//tIg44RcdnD3cT39vZmdDkwExlYLvNDg4Csi0kiN2mpmZkdRusE2OjqIrPpzkiw7mNkPgLMore8+T6m6bUvSdy/WsoOIJBJi2eH4HY9OFHPmLL8n61azV4DPAuuiHz3p7t/p6TlJznYYAvxPYC93/9DMbqNUbLEX6l4sGZkZFwyflPc0pGAuXjar+oMSaOBuh13TPifpssOmwBZm9gmljLeN0j62MdHfq3ux1CzUG00krqnLi919lZn9mlK3ig+Bh939YTNL1L24fAvHZQftwRm7Dwk3e2l6/a6dr09E0oUO1gHMbFvgWGBn4G3gdjObnPQC5RVuZubf/8Ofa5yqtCp9IpJ6afbD1McBr7n7GwBmdifwZRJ2Ly63ZswhmSYrrWfggkXKfKWLYJlvgV9bSYLvcuDAqMrtQ0rHSi4G3idl9+KBCxbVPlNpWcp8pV6aunW8uz9lZrOBZ4D1wL9TWkbYipTdi1/YZUK22UrL2Xvp/YXOTiQfoX4hF3nZQeXFIlJIIfb5jh16RKKYM3/lww3/+NXQCrdLhunoQOnsJytuUuYrXfSGzFflxSLSspp6qxl0W148jZTdi3+yQue2Sle64Sb10sDD1FPLUl4MKbsXrztOW82ks+3vWsT7l1U9Clp6mb7fnxFknFZYdqhUXrxT2ottf5e2mklXod5oInFNHXx7KC/+Mgm6F5eXF/9m7904bdigoP8D0tx2eHAR7047Ie9pSMFsPf3OIOMU+WZu1a1mUXnxHcA3iMqLgdmUGsZV7V4cG6u4/xIiUightpqNHDw6Ucz5v20LC7nVrGJ5sbvfuPEBPXUvLlfk30KSDzPjvX/4Rt7TkILZ6hdVW6Al0uy7HSqWF6t7sYQS6o0mEtfhxT1UMkt58dXqXixZmRnn6zB1ifnHQGc8FznmqLxYRAopxJrvfgO/nCjmPLvm/xRyzTeYOQOU4Uhnx6+dxdVDVHYunZ21KkxBVpHXfJX5ikghhch89xlwYKKY86e1TxYz8zWzcyiVEhtwlbtfZmbbkbJ7cZHXXyQfZsa9/U+u/kDpVSa23xJknCJnvptUe4CZ7UMp8I4E9gMmmtlulM52mO/uuwHzo+9FRAqjwzck+spDksx3T0o96D8AMLOFlLaWHYu6F0sAobIckbgNBf60nST4/gm42My2p7TPdwKlcuLU3YtnzJjBlClTgkxcWoOZaTlKuugN3YsT3XCLWgWdDbwHvEgpCJ/p7v3KHvOWu29bZZzi/kuISKGEuOH2hc//TaKY8+qbzxTzhpu7/xb4LYCZ/SOwkhq6F1+gzfQSc/GyWczrf1Le05CCGdd+W5BxWiHz7e/u7Wa2I/AwcBBwPrDO3aeb2TRgO3c/t8o4xf2XEJFCCZH5Dt/+i4lizrJ1zxUz8wXuiNZ8PwHOdve3zGw6KbsXa21P4rTmK5UEW/Mt8Gsr6bLDVyr8bB2lQ3YS024HqUSvC6mXpj5MPaSFA05s5OWkCYxeO7vQ2YnkQ5lvYKPXzm7k5aRJKPOVemn2fb7dlRdfRMruxUX+LST50JqvVNIb9vkm6V5cXl78MfCgmd0X/XWq7sXKcKQSvS6kXpr6MHW6Ly9ObdGAr9XyNGlhh6y9g+8Mr7pRRnqZf1t2e5BxivypqurBOpTKiw8xs+2jVkITgGHR3001s+fM7Jqo0WYXZjbFzBab2eK7P1gaaNoiItVtcE/0lYcs5cXTUfdiEamTEEUW2261a6KY89Z7rwRZ+zKzHwOXADu4+5s9Pbbm8mJ3X1t2QXUvlpqYGa/uc3je05CC+cKf5gYZp5H7fM1sGHA4paKzqpLudigvLz4BOEjdiyWUUG80kbgGJ3yXAucCv0/y4CzlxTeoe7Fkpa1mUkmoRK1Rux3M7Bhglbs/m3TuWcqLT003PWW+UpleF1IvSW+mlZ87Hpnp7jNjj5kHDKzw9AsoHTR2RJq5NbSB5tvfPEwpjnSy7c2P8p9nac1XOvvc1XOD3HDbfPMdE8Wcjz5aXvO1zGxfSq3UPoh+NBRoA0a6+5pun6fuxSJSRCGC72c3H5Yo5vz1oxXBPn6Z2evAiCC7HUJ5/4qzG3k5aQJ9p17JRy8tyHsaUjCb7zkmyDhFvp+gzFdECilE5rvpZkMSxZz1H69q+I2HhgZf+ZSZTYkv6IvoddF7JCkvlvpQG2epRK+LXkLBV0QkBwq+IiI5UPDNj9b1pBK9LnoJ3XATEcmBMl8RkRwo+IqI5EDBt8HM7Cgze9nMXjGzaXnPR4oh6gbTbmZVj2aV1qDg20Bm1ge4EhgP7AVMMrO98p2VFMS1wFF5T0IaR8G3sUYCr7j7Unf/GLgFODbnOUkBuPsi4C95z0MaR8G3sYYAK8q+Xxn9TER6GQXfxqp0eIf2+on0Qgq+jbUSGFb2/cZDl0Wkl1Hwbayngd3MbGcz2ww4Gbg75zmJSA4UfBvI3dcDU4GHgJeA29z9hXxnJUVgZrOAPwC7m9lKM/u7vOck9aXyYhGRHCjzFRHJgYKviEgOFHxFRHKg4CsikgMFXxGRHCj4iojkQMFXRCQH/x+IsmbhZuWsGQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = sns.heatmap(X_mcar, mask=R_mcar.numpy()==1, linewidths=0.005, linecolor='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Details on all available specifications" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main function `produce_NA` allows generating missing values in various ways. These can be specified through different arguments: \n", "\n", "`produce_NA(X, p_miss, mecha = \"MCAR\", opt = None, p_obs = None, q = None)`\n", "\n", "## Mechanisms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### MCAR\n", "\n", "Missing Completely At Random values are generated using only the desired proportion of missing values `p_miss`, i.e. each value have the same probability `p_miss` of being missing. Therefore, we generate missing values using a Bernoulli distribution of parameter `p_miss`. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percentage of generated missing values: 43.5 %\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAD7CAYAAADNT5fNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAeIElEQVR4nO3dfbQU1Znv8e+jxuBrwCiIgESRhS8kehPHaJwoKL4RQq5ex8ioSUy8mIw6jnlFs+ZqJssMieNE5+pKIMYxowaDKDfEwRdEXkyueEUnGtRhIihyADkSNb5HOee5f3QR+1TXOaeqa3dXdZ/fZ62zpPtU79qyujZP7drPfszdERGR5tqu6A6IiAxEGnxFRAqgwVdEpAAafEVECqDBV0SkABp8RUQKkGvwNbOTzWy1mT1jZjNCdUpEpN1Zvet8zWx74L+AE4AO4BFgmrs/Fa57IiLtaYccnz0CeMbd1wKY2W3AZ4BeB18zU0aHiKTi7pa3jXe3rE015rxvz/1znyurPIPvCGB91esO4OPxg8xsOjAdYNasWUyfPj3HKaXdmBnKspQ4s6aPhU2XZ/BN+tupuYrcfTYwGyqR7/nnn5/jlNKOBsKFJgXp7iq6B73KM/h2AKOqXo8ENvb1AUU4EqfIV5IE+we5a2uYdhogz+D7CDDWzPYDNgBnAn/d1wcU4UgSfS+kUdy7i+5Cr+oefN19q5ldCNwLbA/c6O5P9vWZdR87rt7TSZsa/egDvP3UA0V3Q0pm0MGBxoru8g6+dS81q+tkWu0gIimFWO3wzvrHU405O446tKVWO2TWceTEZp5OWsDIFUt4/Tt9zlbJALTr5T8P01CbPnATESm3Es/55pp2MLPngNeALmCrux/ez/GadhCRVEJMO/xpzYpUY877xxzZktMOE919S5oDtaRI4rTUTJIEWwFT4gduTZ120JIiSaLvhTRMiacd8g6+DtwXTSfMirLZeqhOL77mqAP5wrgROU8p7WTwTYsV+UqNcJFveR+45Z3z3cfdN5rZUGARcJG7L+/jeF1lIpJKkDnfp5ekm/M9aGJrzfm6+8bov51mNp/KTme9Dr6bPnlsntNJGxr+4DJeveCUorshJbP79XeHaajE6cV59vPdBdjO3V+L/rwI+Ad3v6ePzyjyFZFUgkS+T9ybLvL9yEktFfkOA+ZHczM7AD/va+AFrXaQWlrtIElCzfm6l3fON8/eDmuBQ7N8Rk+1JYm+F9IwbbzaIRNFOBKnyFeSaJ1vYIpwJIm+F9IwrRz5mtmNwBSg093HR+/tAfwC+BDwHHCGu7/cX1u/2fu0PH2VNnT0C3eyZfIxRXdDSmbPhb0umsqm690w7TRAmtLxNwEnx96bASx297HA4ui1iEi5dHen+ylAqqVmZvYh4K6qyHc1MMHdN5nZcGCpu49L0Y4m90QklRBLzd5+aE6qMWfQUdNaZqnZMHffBBANwEN7O7A6vfgru/0FJ+58QJ2nlHZ06uY5bP++fYruhpRM17t9loNMbyA/cItXL/7Ra480+pTSYoJdaCJxbTj4bjaz4VXTDp1pPvTk/pPrPJ20q0PWLtRSM6kRLMmixA/c6h18FwCfB2ZG//1lmg8dsnZhnaeTdqalZtIwLb7UbA4wAdjTzDqAy6kMunPN7EvA88BfpTnZmvEn1N9TaUtjVi1S5Cs1lGQBuPu0Xn51fNaTjVm1KOtHZABQ5CsN08qRb0iKcCTOzHjmkBOL7oaUzAFP3hemoVaOfENShCNJgl1oInGtHPn2kl58BfA/gRejwy5z936fpinylTgz482bLyu6G1IyO5/zvTANbQ27mbqZbQ+sBDa4+5Q8baWJfG8CrgP+Lfb+D939n7KcTJGvJAl2oYnEhY98LwaeBnbP21CaB27Lo/Ti3J4//LgQzUgb2XflA6z9yKSiuyEls/8T94dpKOCcr5mNBD4FXAl8NW97aTbW6c2FZvaEmd1oZkN6O8jMppvZSjNbeeuLG3KcTkQkI+9O9VM9TkU/0xNauwb4JhBkRK93Y51hwBYqpeO/Cwx39y+maEeTviKSSoiNdd6aPzPVmLPTqTP6PJeZTQEmu/vfmNkE4OvNmPOt4e6bqzr1E+CulJ+r53TSxsyMm4efVXQ3pGTO2XRrmIbCzfkeDUw1s8nAIGB3M7vF3c+ut8G6Bt9t+zpEL08FVqX8XD2nkzYX7EITiQu02sHdLwUuBaiKfOseeKH+9OIJZnYYlWmH54Dz05xMG+tI3CFrF3L/0DOK7oaUzKTOuWEaKvHddqo532An05yviKQUZM53zuXp5nynfadlNlOvi+Z8Jc7MeOOqLxXdDSmZXb7x0zANKb24QnO+kiTYhSYS1+LpxaOoZLftTWV922x3v7aeCsYPDT81b3+lzRy1ab7uiKRGsECtqytMOw2QJvLdCnzN3R8zs92AR81sEfAFKhWMZ5rZDCoVjL/VV0NHbZqft7/ShnRHJA3TytMO0ZKybcUyXzOzp4ERwGeorIIA+BmwlH4GX0U4Emdm/PMorfOVnr66PtDywxIPvpnSi6NMt/8GPEysgjGQWMG4Om1v9uzZ+XorIpJFyvTiIqReamZmuwLLgCvd/U4ze8XdB1f9/mV373WPh+gYhb4ikkqIpWZvzr4k1Ziz8/QflnOpmZm9D7gDuNXd74zezlzBWNMOEmdmHD78k0V3Q0pm5aYHwzRU4mmHNKsdDPgp8LS7/3PVrzJXMNaDFUkS7EITiWvx1Q5HA+cAvzOz30bvXUYdFYyXDD293n5Km5rYOY+7hp5ZdDekZKZ03hamoRJHvkovFpFSCjLne+2X0835Xvzjcs75hvK7/T7VzNNJC/jws//O9NG6I5KeZq+bF6ahEj9naurg++Fn/72Zp5MWEexCE4kr8bRDnvTiK8hYwVirHSTOzDh1308X3Q0pmfnP/ypMQ93lHXPypBdDxgrGWu0gSYJdaCJxrbzaoY/04swU+UqcmfH9fZVeLD196/kw6cXeytMO1WLpxUdTqWD8OWAllei4ZlezqApoUiVQESDchSZSo8WnHYA/pxffAfydu79qZj+iUrl4WwXjq4GaCsbuPhvYtqlDef8mpBBmpjJCUiNcGaEWj3yT0ovrrWAsItI0rRz59pZeXE8FYz1wkyTBohyRuK0t/MCN3tOLp2WtYKwHbhJnZjzxoSlFd0NK5iPPBbqRbuVpB3f/NZAUsva5pjeJIl9JEuxCE4lr5WmHkBT5SpyZsfrAk4ruhpTMuP+8N0g7Lb3UzMwGAcuB90fHz3P3y+spoKnIV5KEutBEarR45Psn4Dh3fz1a9fBrM7sbOI2MBTQXD+1310kZYI7vvF13RFIjWKDWyoOvV66M16OX74t+nDoKaB7feXud3ZR2pjsiaZhWTi8GMLPtgUeBA4Dr3f1hM+tRQNPMEgtoVlPkK3GKfCVJqH+QvZUjXwB37wIOM7PBwHwzG5/2BEovlv4o8pWGafXBdxt3f8XMlgInk7KAZnV68Q9Gn13evwkpxLeev1WRr9QIN+fb2qsd9gLejQbenYBJwPepo4CmNlCRJIp8pWFaPPIdDvwsmvfdDpjr7neZ2UNkLKCpCEfizIyOIycW3Q0pmZErloRpqJUHX3d/gso2kvH3/wAc34hOiYiE4F0tPO0Qkm4vJUmwKEckrpUj35CWDVOVWunp2M3zNB0lNbTUjD7Ti68gYwHNYzerSq3U0h2RNEwrD770nl4MGQtoKsKRODPjqlGq4SY9fWN9oJVRgaZ8e6vinqfNPOnFmSnCkSTBLjSRGN8a7IFbYhV3d3+q3gbzpBefQooCmtU2Hj2h3n5Km9rnN0t59XxtKSk97T4r0E53gcbePqq41z34WpapgG3pxcBFVOZ6t/BeAc3h7l5TQDOWXvyxejsqIgOLu+e+VX75ryakGuCG3L409bmiKu7LgfHu/mp9PcuRXlw919tXAU1VL5a+mBlv3vC1orshJbPzeVeHaShl5JuwB83saOyKH9ejinuertWdXqwCmhJKsAtNJCbtUrNYkJgoqYp7HnnSi29WAU3Jy8z0vZAa4TbWCdNMb1Xc88iTXnxO1pMp8pUk+l5Io/jWYE0lVnHvL7ehL03NcPs/w6Y183TSAv775jlcP1LrfKWnCzrCLD8MVTm+jyrudWvq4Csi0lTl3Vcn/VKzaM53JbDB3afUWb1Yk3sikkqIpWYvnnBsqjFnr0XLmj73lSXyvRh4Gtg9ej2DjNWL9WBF4syMJUO14ZL0NLEzzD4woaYdGiFthttI4FPAlcBXo7czVy/WgxVJEupCE4nzrvKOOWkj32uAbwK7Vb2Xqnpx9eLl7ww7hM8OHpWju9JuDlx9DxeMPqPobkjJXL9ubpB2WjryNbMpQKe7P2pmE7KeoHrxspn55ZufzNxJaW+hLjSROO9u7cj3aGCqmU0GBgG7m9ktpKxeXG3BsDPz9VbaztTNt+lZgNQIt5l6kGYaIuvGOhOAr0erHa4C/lD1wG0Pd/9mP5/XVSYiqYRY7bDhqONSjTkjHnqg1Ksd4mai6sWSk9KLJclAiHyz7mq2lMqqhrqqF2u1gyTR90IapbsNVjsEsemTxzbzdNIChj+4TJGv1AgX+WrwFRFpurYYfBPSi68gY/Xi4Q8uq7ef0sY07SCNUuabqjzpxZCxevGa8SdkOJ0MBGNWLeKS0Z8tuhtSMj9c94sg7bR85NtLenFmY1Ytqvej0sZCXWgicQFWqzVMnvRiSFG9uDq9eNasWUyfPj1+iAxgWmomSUJNRXWVeLVDv0kWUXrxZHf/m1iSxTBSVC+OtaWrTERSCZFksfrAU1KNOeP+8+5SJlkkphe7+9nbDuirenG1P0w9pu6OSnv64ILlehYgNUJNUZZ5zjdPevGfqxeb2SXAx929z80bFPmKSFohIt+nx05ONeYc9PuFpYx8e/MDVS+WvMyMjUdPKLobUjL7/GZpkHbaJvLNfTJFviKSUojI93f7fTrVmPPhZ3/VUpFvZq9+5eRmnk5awO4/uocVw08tuhtSMkdumh+knTLfbCvyFZFSChH5/nb01FRjzmHrFpQz8jWz54DXgC5gq7sfXk/14ucPPy5PX6UN7bvyAT0LkBrBNtYpcZLFdhmOnejuh7n74dHrbdWLxwKLo9ciIqXhnu6nCHnmfDNXL9535QM5TiftShvrSKN0lzjyTTv4OnBfNGc7KyqKmbl6sdKLJU7pxZIkWHpxd5ab++ZK9cDNzPZx943RALsIuAhY4O6Dq4552d2H9NOOrjIRSSXEA7cV+5yWasw5cuOd5Xzg5u4bo/92mtl84AjqqF6sCEfizIxzR/+PorshJfOv6+4I0k5LTzuY2S7Adu7+WvTnE4F/ABYAn6dSSPPzwC9TtJWvt9KWQl1oInFlXu2QJvIdBsyPBs4dgJ+7+z1m9ggZqxdvmayNdaSnPRcu1x2R1AgVqJW4eHH/g6+7rwUOTXg/c/XiPRcuz3K4DBC6I5JGccr73WpqerEiHIkzM5YMPb3obkjJTOycF6SdrS0+7RCMIhxJEupCE4lr+ci3l/TiK8hYvfiN6y6ov6fSlna58HrdEUkNzfn2NNHdt8Tey1S9eJcLr89wOhkodEckjdLykW8or/+9SoRLT7t+9xe89etbiu6GlMxOf3l2/welEDLyNbOTgWuB7YEb3H1mnvbS5t5tSy9+NEoX3uZCM3vCzG40s8TsNjObbmYrzWzljY+tydNXEZFMurBUP/0xs+2B64FTgIOBaWZ2cJ6+5UkvXo2qF4tIg4RIL/7V3tNSjTmffmFOn+cys6OAK9z9pOj1pVEf/7HevtWdXuzuf160m7Z6sR6sSJw21pEk4R64pWunegOwyOxoA7FtRgDrq153AB/P07e604urqxcDpwKrUrSVp6/SpvS9kEZJ+896NNDO7uOQpC9prqghT3rxzVmrF284akL9PZW2NOKhpdw/9IyiuyElM6lzbpB2Aj5w6wBGVb0eCWzM06BquIlIKYWY8503/KxUY87pm27tb853B+C/qGypsAF4BPhrd3+y3r41dalZx5ETm3k6aQEjVyzRnK/UCLaZepBWwN23mtmFwL1UlprdmGfghSYPviNXLGnm6aRFaM5XGqU74FcryuDtM4s3i7TpxYOBG4DxVOZ4v0hlqVmm6sUvnzmh/p5KWxpy21JeOu3YorshJbPHncuCtJN2tUMR0q7z/RnwoLvfYGY7AjsDlwEvuftMM5sBDHH3Pgtoas5XRNIKMed7yz5npxpzzt54S/nKCJnZ7sAxwBcA3P0d4B0zy1y9WHN7EmdmvHS6Il/paY95gSLf8ga+qaYd9qeyc9m/mtmhwKPAxdRRvVgkSagLTSSu1Xc12wH4KHCRuz9sZtcCM9KeILZ4WaGv9GBmPHvopKK7ISWz3+P3B2mnq8SRb5qNdTqADnd/OHo9j8pgvDmqWkza6sUiIs3UnfKnCGlquL1gZuvNbJy7r6ayyPip6EfViyW3UFGOSFyrTztAZRezW6OVDmuBc6lEzZmqF+uBm8RpYx1JEipQK3EJt9S7mv0WODzhV5mqFyvylST6XkijtEPkG4QiHIkzM24eflbR3ZCSOWfTrUHaCZVe3AiqXiyFC3WhicS1+jrf3tKLTyJj9eLOScfU31NpS0PvX647Iqmh6sXvuRa4x91Pr0ovPomM1YuH3r+8/4NkwNEdkTRKSw++faQXZz7ZUweckvkz0t4OfuZu/mWk5nylp7/tCDMVVeZ7qjzpxVCpXvw5YCXwtaRdzZReLP0JdaGJxJV5zrffXc3M7HBgBXB0VXrxq8B1ZKxeTLn/IZICaJ2vJIm+F7mHzn8cnW5Xs0vXNX9Xs7rTi919s7t3uXs38BPgiEZ1UkSkHt14qp8i1J1erOrFEoq+F9IoLf3ALZKUXvwvWasXP3PIiXV2U9rVAU/ex4376IGb9PTFje3/wE3Vi0WklELM+V4xOl314ivW9V29uBGUXiyF0gM3SRJqKmprieO9NOt8x1EplLnN/sD/Av6NjAU0NbcnSfS9kEYp79Cb7oHbauAwADPbHtgAzKdSzWJxVQHNGfRTw23J0NNzd1jay8TOeYp8pYbSi2sdD6xx93X1FNCc2Dkva/9kAFDkK41S1DKyNLIOvmcCc6I/pyqgWU0RjsRpzleSBNtMPUgrjZF68I2WmU0FLs1yAqUXS38U+UqjtMu0wynAY+6+OXq9eVuiRV8FNKurF2+eMKHM/xBJAfZetow1408ouhtSMmNWLQrSTleJY98sg+803ptyAFhAxgKaey9blqlzMjCEutBE4lo+8jWznYET6JnFNhMV0JSczIx3XlxTdDekZHbca0yQdrzVI193fxP4YOy9P5CxgKaISDO1fOQbih6sSJJQUY5IXDstNctl3ceOa+bppAWMfvQBTUdJDS01o8/04sFkLKA5+tEH6uymtDPdEUmjbC3x8JsnvfhcMhbQVIQjcWbGRaPPKLobUjL/e93cIO20/AO3KtXpxZlPpghHkoS60ETi2umBW3V6MaQooFntj+dpMb309IEbFumOSGqEm/Mt73crT3rxj6gUztxWQPNqoKaAptKLpT+6I5JGaVbka2ZXAZ8G3gHWAOe6+yt9fabu9OKqNGPM7CfAXUkfqk4vptwPH6UAZsbNw1VGSHo6Z1OYMkJdzburWgRc6u5bzez7VILUPnd5rDu9WAU0JZRQF5pIXLPW+br7fVUvVwD9bl6eJ734B1kLaGpuT+K0paQkafacb8L06Ozorr0eX6Tn8txEedKLz8naI0W+kkTfC2mUtHO+senRRGZ2P7B3wq++7e6/jI75NrAV6Pd2rqkZbo+PntLM00kLOHTdXYp8pUa4MkLhvlvuPqmv35vZ54EpwPGe4kvd1MFXRKSZmrXUzMxOpvKA7dhopqBfaed8LwHOozK/+zsq2W07k7F68aHrEhdEyACnaQdplCaudrgOeD+wKPo+r3D3L/f1gTR7O4wA/hY42N3fMrO5VJItDiZj9WLdXkqcmfHt0dOK7oaUzJXr5vR/UApNXO1wQNbPpJ122AHYyczepRLxbqSyjm1C9PtU1YsV4UiSUBeaSFxLpxe7+wYz+ycq1SreAu5z9/vMLFX14uolHNccdSBfGDciXO+l5Q2+aTE3jFCShfR03oYwa7/LnF5s/U0FmNkQ4A7gs8ArwO3APOA6dx9cddzL7j6kn7bK+zchIqXi7rlvlSfvOznVmLPw+YVNvy1PM+0wCXjW3V8EMLM7gU+QsnpxNc35SpySLCRJsCSLEn+30gy+zwNHRllub1HZVnIl8AYZqxdrzleS6HshjdLSpePd/WEzmwc8RiVz4z+oZILsiqoXS05mxpdH9/vVkQHmx+tuD9JOmWu49TvnG/RkmvMVkZRCzPkeP/LEVGPO4o77SjnnG8xVo/RUW3r6xvpbdUckNcqYXhya0otFpG2VealZnvTiGWSsXvyN9dq3VWrpgZs0ShPTizPLk14Mql4sOZkZb1zT71bQMsDs8nezgrTTDtMOSenFH8p6MkU4kiTUhSYS19KDbx/pxZ8gRfXi6vTiqw8Zy+dGDQ/6PyCtba97lvPajNOK7oaUzG4z7wzSTpnvtvOkFy8CtvBe9eLh7l5TvTjWVnn/JkSkVEIsNTtin2NTjTn/b+OyUi41S0wvdvdbth3QV/Xiak+PPaXefkqbOuj3d5c6OpFiNLuGWxHqTi+up3rxQb+/u+6OSvvSswBplC4v76aSedKLb8havfjJ/Sfn6qy0n0PWLlTkKzUGwsY6Si8WkVIKMed76N6fSDXmPP7C/y3lnG8w84epXIz0dOrmOdpMXWpoM/XQJ1PkKyIphYh8xw87MtWYs2rzinJGvmZ2MZVUYgN+4u7XmNkeZKxe/Nioqbk6K+3no+sXlHpeTooxEFY7bNffAWY2nsrAewRwKDDFzMZS2dthsbuPBRZHr0VESqPLu1P9FCFN5HsQlRr0bwKY2TIqS8s+Q8bqxR9dv6Defkob01IzaZTuEt9VpRl8VwFXmtkHqazznUwlnThz9eJzP3AEx+08NkjHpT2cs0n7+UqtgTDtkOqBW1Qq6ALgdeApKoPwuapeLCKNEuKB25g9P5pqzFmz5bFyPnBz958CPwUws+8BHah6sQRgZtw/9IyiuyElM6lzbpB2yhz5pl3tMNTdO81sX+A04ChgP1S9WAIIdaGJxHV5V9Fd6FXaJIs7ojnfd4EL3P1lM5tJxurFdw09s79DZICZ0nkbr31NSxClp92uDvNwvsx320qyEJFSCjHnO3KP8anGnI6XVpVzzjeUMv8rJMUwMy4brbRz6el76+YEaafMY44iXxEppRCR7/DBB6cacza98lQ5I99e0ouvIGP14hdPPiZHV6Ud7XXPcq4apY11pKdQlc7LvNohTRmh8cBtVNKL3wHuAb4CnAW8nqV6sSJfEUkrROS71wfGpRpzXvzj6lJGvr2lF2dW5vkXKYaZ8eXR/S6UkQHmx+tuD9JOmcecfjfWoZJefIyZfTAqJTQZGBX97kIze8LMbowKbdYws+lmttLMVs6ePTtQt0VE+tftnuqnCHnSi2ei6sUi0iAhph2G7HpAqjHn5defCTLtYGZfB64C9nL3LX0dW3d6sbtvrjphqurFC4YpyUJ6mrr5NtaMP6HobkjJjFm1KEg73U184GZmo4ATqCSd9X98ysi3Or34PirpxYO27WpmZpcAH3f3PkdXRb4iklaIyHf3XfZPNea8+sba3OeKCg1/l8pWC4cHiXxJTi++OWv1YkU4Ejdm1SJePmNC0d2Qkhkyd2mQdpq1UbqZTQU2uPvjafewUZKFiJRSiMh3p51Gpxpz3n77+fOJ9h2PzHb3HisEzOx+YO+Ej38buAw40d3/aGbPkSLyberg+8pZx2nwlR6G/HwJfzxPd0TS0wduWBRk8B00aN+0g2/d5zKzD1MppfZm9NZIYCNwhLu/0OvnFPmKSBmFGHzfP2hUqjHnT2+vD5ZkkTby1cY6Uigz4+2nlxbdDSmZQQdNCNJOmcccRb4iUkohIt8ddhyRaszZ+s6GpqcXN3XwlfeY2fT4hL6IvhcDR5r0YmmM6f0fIgOQvhcDhAZfEZECaPAVESmABt/iaF5Pkuh7MUDogZuISAEU+YqIFECDr4hIATT4NpmZnWxmq83sGTObUXR/pByiajCdZraq6L5Ic2jwbSIz2x64HjgFOBiYZmYHF9srKYmbgJOL7oQ0jwbf5joCeMbd17r7O1SqQn+m4D5JCbj7cuClovshzaPBt7lGAOurXndE74nIAKPBt7mSNu/QWj+RAUiDb3N1AKOqXm/bdFlEBhgNvs31CDDWzPYzsx2BM4EFBfdJRAqgwbeJ3H0rcCFwL/A0MNfdnyy2V1IGZjYHeAgYZ2YdZvalovskjaX0YhGRAijyFREpgAZfEZECaPAVESmABl8RkQJo8BURKYAGXxGRAmjwFREpwP8HP41PyiayilsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Sample mcar data -----------------------------------------\n", "X_miss_mcar = produce_NA(X_complete_cont, p_miss=0.4, mecha=\"MCAR\")\n", "\n", "X_mcar = X_miss_mcar['X_incomp']\n", "R_mcar = X_miss_mcar['mask']\n", "\n", "print(\"Percentage of generated missing values: \", (R_mcar.sum()).numpy()/np.prod(R_mcar.size())*100, \" %\")\n", "\n", "ax = sns.heatmap(X_mcar, mask=R_mcar.numpy()==1, linewidths=0.005, linecolor='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### MAR\n", "\n", "Missing At Random mechanism values are generated by using a logistic model. A subset of fully observed variables (with *no* missing values) is randomly selected. The remaining variables have missing values according to a logistic model (depending on the fully observed variables only) with random weights, re-scaled so as to attain the desired proportion of missing values on those variables.\n", "\n", "There is an additional argument: \n", "\n", "* `p_obs`: the proportion of fully observed variables that will be used for the logistic masking model. Note that at least one variable is chosen to be observed\n", "\n", "Note that there exists other ways to generate missing values as described in the Rmarkdown [How generate missing values?](https://rmisstastic.netlify.app/how-to/generate/misssimul) by using patterns as in the first definition of Rubin (1976)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percentage of generated missing values: 21.0 %\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAD7CAYAAADNT5fNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAf6UlEQVR4nO3deZhU1Zk/8O9XTVxxQGVfWiMMasy44A9NNIqibCLuCooxan5k0zEmDiHGmTgZzeAyUWdMMnRc4gIigiiiCIgC0YkLmHFFRlwamq1RcQQ1Kt3v/FG3Y9Xt6q5zq07XPVX9/TxPP1rdVecefbgvp9973vPSzCAiIuW1XdoTEBHpiBR8RURSoOArIpICBV8RkRQo+IqIpEDBV0QkBSUFX5IjSK4kuYrkJF+TEhGpdix2ny/J7QH8D4ATANQDeB7AODN7zd/0RESq0w4lfHYwgFVm9hYAkJwO4GQArQZfkqroEBEnZsZSx/j83becYs6X9vpKyddKqpTg2xvAmqzX9QAOj7+J5AQAEwDgxsED8e3+vUq4pFSbLtOexMdzb0x7GhKYXUZflvYU2l0pwTff3xQt/pYxs1oAtUBm5XvZcytLuKRUo45wo0lKmhrTnkGrSgm+9QD6Zr3uA2BdWx/QORISRxIfjD8u7WlIYDrf84SfgRq3+RmnHZTywG0HZB64DQWwFpkHbueY2attfEbRV0Sc+Mj5frbuVaeY8+VeX62cnK+ZbSN5MYD5ALYHcHtbgTf6TLGXkypFEp88OyPtaUhgdj78LD8DNTX5GacdFL3yLepiWvmKiCMvK981L7qtfPseVDkr32Jo5StxJLHl8lPSnoYEptMND/oZqEofuImIhM2qNO1A8h0AWwA0AthmZocVeL+WviLixEfa4dM3n3GKOTvue0RFph2ONbN3Xd64cr/hHi4n1WTg6/PxdI/T0p6GBObIDQ/4GSjgB25lTTsMfH1+OS8nFcLbjSYSF3DaodTgawAWROmEKVE1W47s8uIpU6ZgwoQJJV5SqglJ/LnfSWlPQwJzyOqH/QwU8AO3UnO+vcxsHcluABYCuMTMlrbxfuV8RcSJl5zviifdcr77H1tZOV8zWxf9s4HkbGROOms1+NYNUhmp5KpZ/gQ++PbQtKchgen8h0V+Bgq4vLjo4EtyVwDbmdmW6N+HAfhlW5+pWe6pXluqircbTSSuSh+4dQcwm2TzONPM7LG2PqAiC4kjifXfPCbtaUhgev5xiZdxzMLN+ZZytsNbAA5K8pkoUIvk8HWjibRQxbsdEtHKV+JI4uN7rkx7GhKYXcZf7WegKk07JKaVr+Tj7UYTiavklS/J2wGMBtBgZgdG39sDwH0A9gbwDoCzzGxzobGe7HZGKXOVKnRsw0xsOEY5X8nVY4mnVFTj537GaQcureP/AGBE7HuTACwyswEAFkWvRUTC0tTk9pUCpyILknsDmJu18l0JYIiZrSfZE8BiMxvoMI6SviLixEeRxV/+dK9TzNnp6+Mqpsiiu5mtB4AoAHdr7Y3Z5cXf2X0wjt+lf5GXlGo0dsM0bP8ldbSWXI2ft9kO0l1HfuAW715864fPtfclpcJ4u9FE4qow+G4k2TMr7dDg8iFtNZM4ktjykzFpT0MC0+nf5ngZxwJ+4FZs8J0D4HwAk6N/PuTyIW01k3x83WgiLVT4VrN7AQwBsBfJegC/QCboziB5EYDVAM50uZhWvhJHEp9tejPtaUhgvtx1Xz8DVXLawczGtfKjxEdRaeUr+Xi70UTiKnnl69PyPsrtSa5B9XPwWv+RaU9DAnPAqnl+Bqrkla9Pg+qV25OWvN1oInGVvPJtpbz4KgD/H8Cm6G1XmNmjhcbaeqVTalg6kN2uvh8f/fbitKchgdn1B7f4GWhbuIepF6xwI3k0gK0A7ooF361mdkOii6nCTUQc+ahw+2Tur51izs6jfxxehZuZLY3Ki0um3Q4SRxKvD4wfHSId3X4r2+zL4C7gnK/LwTqtuZjkSyRvJ9mltTeRnEByGclltbUtmhuLiLQfa3L7SkGxB+t0B/AuMq3j/wVATzO70GEcLX1FxImXtMPsyW5ph1MnhZd2yMfMNjb/O8nfA5jr+LliLidVjCRu73Vu2tOQwFy4bqqfgSp5t0M+zec6RC9PBfCK4+eKuZxUOW83mkic590OJLcHsAzAWjMbXcpYxZYXDyF5MDJph3cAfNflYi/0VZGF5Dp0zRw80u3stKchgTmx4T4/A/n/bftSACsA7F7qQMWWF99WzMUOXaMiC2nJ240mEudxtwPJPgBOBHANgB+XOp66F0uqSGLrL8amPQ0JzG7/PN3PQI7BN7vpQ6Q2Oos8200AJgLo5GNq6l4sqfN2o4nEOT5wy276kA/J5irf5SSH+JiaS863L4C7APQA0ITM3wg3F9PBWN2LJe7YhpmYp5yvxIz0lYpqbPQzDnAkgDEkRwHYCcDuJO8xs/HFDuhSXtwTmX28L5DsBGA5gFMAfBvA+2Y2meQkAF3M7KcFxlLeQUSceNnne8dEt32+F1znfK1o5Xt5u+92iLaUNTfL3EJyBYDeAE5GZhcEANwJYDGANoOvcr4SRxLX9tM+X8n109Weth9WS3lxVOl2CIBnEetgDCBvB2OVF4tIatqhvNjMFpe66gUcy4sBgORuAJYAuMbMHiD5gZl1zvr5ZjNr9YyH6D1a+oqIEx9ph49rL3OKObtMuDHM8mKSXwIwC8BUM3sg+nbiDsZ1g44rfqZSlWqWP4FBPY5KexoSmOUbnvIzUMBpB5fdDkSmqGKFmf0660eJOxjXLH+iyGlKNfN2o4nE+dvt4J3LyvdIAOcBeJnkf0ffuwJFdDCery1FEjO84T7M6t5aj1bpqE7feK+fgQJe+TrnfL1cTDlfEXHkJed78/fccr6X/meYOV9f1L1Y4gbVz8EFNaenPQ0JzB11s/wMFPD2VnUvltR5u9FE4gJOO5RSXnwVEnYwvrunNtNLrvPWT8WYviemPQ0JzJw1j/gZqCnclW8p5cVnIWEHY+V8RcSVl5zvtRe45Xx/ekd4Od82yosTU3mxxJHE1f3OSXsaEpgrV0/zMo5VctohW6y8+EhkOhh/C5m2Gj/Jd6pZnnMyRXL4utFEWgg47eAcfKPy4lkAfmRmH5L8HTKdi5s7GP8bgBYdjLPPybyj9/hw/09IKi5cN1VthKQFf22EKnzlm6+8uNgOxiIiZVPJK9/WyouL6WCsLrWSj3q4SbvZVp3lxeOSdjB+Y/9hRU5TqtWAFQvwfO+T056GBOb/rS14VIybSk47mNlTAPJtw2hzT28+A1YsSPoR6QC83WgicZWcdvBJW80kjiRe3kdFFpLra2/7KbKo6K1mJHcCsBTAjtH7Z5rZL4ppoKnuxZKPrxtNpIUKX/l+CuA4M9sa7Xp4iuQ8AKcBWJTVQHMSCvRwU5daiRvZcB/u7KWyc8l1vq+H85UcfC2TK9gavfxS9GUoooGmt3bQUlW83WgicRV+mDpIbo/MmQ79AfzGzJ4lmdNAk2TeBprZtPKVuJEN9+F2rXwlxte2VKvklS8AmFkjgINJdgYwm+SBrhdQebEUov3f0m4qPfg2M7MPSC4GMAKODTSzy4uvqTk33P8TkoorV09Tzlda8JfzDXe3g8uRkl0BfB4F3p0BLABwLYBjALyX9cBtDzObWGAsBV8RceLjSMktPxjpFHM6/XZeeEdKAugJ4M4o77sdgBlmNpfkn5Cwgab2+UocSbx90PFpT0MCs8+Lj/sZqJLTDmb2EjLHSMa//x6Aoe0xKRERH6wx3LRDWSvcVGQh+Xhb5YjEVfLK16cF3c4q5+WkAgxrmIFhfUakPQ0JzIL6x7yMU9FbzdooL74KCRtoDmuYUdpspSr5utFEWqjk4IvWy4sB4MYkDTQXdSv4TE46mKEN9+NX/bTVTHJdsdrXVjM/w7SHUsqLExvacH8xH5Mq5+1GE4mxbeFG31LKi0fCoYFmNm01kziS+GD8cWlPQwLT+Z4n/AwUbuwtXGSR8+aovBjAJcjket/FFw00e5pZiwaasfLiQaVOWEQ6Bh9FFpvPHOIU4LrcvzjIIou/yi4vzs71ttVAM7u8GEWmK6R6kcRHNxXsQCUdzK4/muJnoIBXvi67HeLlxccDuLaYBpra5yv5eLvRRGJ8bTUj2RfAXQB6IBPSa83s5lLGLKW8+O6kDTSf7nFaKXOVKnTkhgewadjRaU9DAtN1wVI/A/lb+W5D5rnWCyQ7AVhOcqGZvVbsgIlyvqXSwToi4spHzve9E49xijl7PrIk0bVIPgTgFjNbWNTEUOYKtxk9zinn5aQCnLVhGm7qq32+kutHa3wdpu5lmBwk90bmvJtnSxlnOx+TEREJUpPbF8kJJJdlfeVtAEFyNwCzAPzIzD4sZWrOaYco57sMwFozG11k92KlHUTEiY+0w6YT3NIOXRcWTjtEFb5zAcw3s1+XOrckaYdLAawAsHv0ehISdi9e3mdMUZOU6jWofg7mq7efxAz31GzXV9qBma1atwFY4SPwAo4rX5J9kOlQfA2AH0cr35UAhmS1EVpsZgMLjKOVr4g48bHy3TjErcii++K2iyxIHgXgjwBexhd7KAoeJtYW15XvTQAmAuiU9T2n7sXZFW7/tOfXcEanmmLnKlXo796Zi+/V6MAlyfWfdX7OgfG18jWzpwB4LVRwKbIYDaDBzJaTHJL0AtkVbiTtl++9nHiSUt183WgicdYUbmGXy8r3SABjSI4CsBOA3UneA8fuxdl0sI7EkcTms4akPQ0JTJcZi72M0x5bzXxJerDOEACXRznf66HuxSLSTnzkfNd+/TinmNP7T0+EfbBOzGQk7F780XUXlHA5qUa7TrwDH92Yd0uldGC7XlZb+E0OqmblW/LFtPIVEUc+Vr6rDxvqFHP6LVtUUSvfxJTzlTiS+o1IWth14h1exqn0B24iIhWpKoJvnvLiq5Cwe7HO85V8fK1yROJC/mW7lPJiIGH34hUDRia4nHQE+78xD39fo/JiyfXvdb7Ki8Nd8Lk20OwD4ERE5cXFXmz/N+YVfpN0OL5uNJE4D8/s2k0p5cWAQ/fi7PLi6/b9W4zv0auE6Uq16fX0YtzTU+f5Sq7x6/2c59vYGG7wLbjVLCovHmVmP4gVWXSHQ/fi2FgBZ2BEJCQ+tpqt3G+kU8wZ+Pq8ILea5S0vNrPxzW9oq3txto3HqVeX5Or+xFI9C5AWfKUoQ875llJe/NfuxSQvA3C4mY0t8HmtfEXEiY+V74oBo5xizv5vPBrkyrc11yXtXqwiC4kjibpBx6U9DQlMzfInvIxTNSvfki+mla+IOPKx8n15n5OcYs7X3n64ola+iX3wLa1wJFfnu57A0u6npz0NCczRG2d5GSfkX7a18hWRIPlY+f53zRinmHNw3ZwwV74k3wGwBUAjgG1mdlgx3YvfPPCEUuYqVWjfVxbitf7a7SC5DljlabdDwEUW2yV477FmdrCZHRa9bu5ePADAoui1iEgwzNy+0uDavfgdAIeZ2btZ31P3YhFpNz7SDsv6nOIUcw6rfzDMtAMy28kWRMFzStQUM3H34ilTpmDCBHUtkC+QxCfLH0p7GhKYnQed7GWcxqYkv9yXl+vKt5eZrYsC7EIAlwCYY2ads96z2cy6FBhHK18RceJj5ftMr9OcYs4R6x4Ic+VrZuuifzaQnA1gMIroXvxizeiSJivV56C6ufhWzWlpT0MCc1fdA17GaQr4gZvLwTq7AtjOzLZE/74QwC8BDIW6F4tIO/Gx8n26xxlOMefIDTODXPl2BzA76kKxA4BpZvYYyeeRsHuxyosljiTG9zs17WlIYO5ZPdvLOAE3L1aRhYiEycfKd2mPM51iztEb7g9y5evNg93HlfNyUgFO2Xgv5ndTGyHJNbzBT3eTbZWc8/V6Ma18RcSRj5Xvou5nO8WcoRvvC3Pl20p58VVI2L1YOV+JI4mPH7057WlIYHYZdamXcULO+SZJOxybXeEWSdS9WK3jJR9fN5pInCHcmFPWnK9WvhJHEh8/9h9pT0MCs8uIS7yME/LK17X2rrm8eHlULtzsYpIvkbydZN7qNpITSC4juay2trbkCYuIuGoEnb7SUEp58Uqoe7GItBMfD9we7jHOKeactOHeMB+45SsvNrOlzT937V68aYS6F0uuro8txcSaNvuuSgd0Xd10L+M0eVzVkhwB4GYA2wO41cwmlzJeweCbp7x4GIBfZncvBnAqgFcKjdX1saWF3iIdkK8bTSTO16/aJLcH8BsAJwCoB/A8yTlm9lqxY5ZSXnx30u7Fbx90fLHzlCq1z4uP4xEVWUjMiZ6KLDw+cBsMYJWZvQUAJKcDOBlA+wXf6GIH5fn+eUkvts+Ljyf9iHQAvm40kbgmf9tbewNYk/W6HsDhpQyorWaSKpLYcvkpaU9DAtPphge9jNPo+L7spg+R2qhpxF/fkudjJQW0sgZfFVlIPr5uNJG4JseQEwXatvbC1gPom/W6D4B1RU8M7uXFnQHcCuBAZKL9hchsNUvUvVgrX4kjiU3DtAtGcnVd4OfhvMfdDs8DGEByHwBrAYwFcE4pA7qufG8G8JiZnUHyywB2AXAFMt2Lmw9TnwTgp20NopWv5OPrRhOJ87XcM7NtJC8GMB+ZrWa3m9mrpYzpstVsdwBHA/h2NInPAHxG8mQAQ6K33QlgMQoEXx0dKHHDG+7T/m9pwde2VNe0g4vo4LA2Dw9LwmXl+xVkTi67g+RBAJYDuBRFdC8WyUf7v6W9hHy2g0vw3QHAoQAuMbNnSd6MTIrBSXYie/OZQ5T0lRx7zFyClfsNT3saEpiBr8/3Mk5jwJlOl4N16gHUm9mz0euZyATjjVHXYrh2LxYRKacmx680uBRZbCC5huRAM1uJTNfi16Kv8wFMjv75UKGx9pi5pMTpSjXytcoRiav0tAOQOcVsarTT4S0AFyCzak7UvXhRt4JvkQ5maMP9qO19btrTkMBMWDvVyzgBt3BTDzcRCZOPIyV/23e8U8z5wZp7wjxS0hcVWUgcSdyqla/EfMfTyte1vDgNKi+W1Pm60UTifO7z9a2U8uLhUPdiKRFJvH/GMWlPQwLj6+F8NTxwy1dePBzqXiweaBeMtJeKDr5tlBcnvtif+52U+DNS3Q5Z/TBu6Kucr+S6fI2n3Q5eRmkfpZQXA5nuxd8CsAzAT/KdaqbyYinE140mEhdyzrfgVjOShwF4BsCRWeXFHwK4BQm7Fz/d44yQ/yKSFBy1cZb2+UoLE9ZO9bLV7F9r3Laa/ayu/FvNii4vNrONZtZoZk0Afo9MjyMRkWA0wZy+0lB0eXEx3YuP2jirtNlKVfJVzSQSV9EP3CL5yov/PWn34tf6jyxymlKtDlg1D7/ro7SD5Pp+ffU/cFN5sYgEyUfO96qac51izlV1U1VeLB0LSbx/moosJNceD/jZ+70t4PWeyz7fgcg0ymz2FQD/BOAuJGygqSILycfXjSYSF27odXvgthLAwQBAcntkOnfORqabRaIGmlr5ShxJvDtKPdwk116P+upeHK6kaYehAN40s7piGmhq5Sv5+LrRROLS2kbmImnwHQvg3ujfnRpoZnu821kJLyfV7viGGThVZecSM3v1w17GCTf0Jgi+0TazMQB+luQCKi+WQnzdaCJx1ZJ2GAngBTPbGL3e2Fxo0VYDzezuxfWHHxfyX0SSgr7PPYkVA7T/W3Lt/8Y8L+M0Brz2TRJ8x+GLlAMAzEHCBpp9n3sy0eSkY/B1o4nEVfzKl+QuAE5AbhXbZCRsoKndDhJHEp+++WzhN0qHsuO+h3sZxyp95WtmHwPYM/a995DZ/SAiEqSKX/n6oq1mko+vVY5IXDVtNSvJqq8OK+flpAL0f3UBrqgZl/Y0JDC/qru38JschBt6Sysv7oyEDTT7v7qgyGlKNfN1o4nEbQs4/JZSXnwBEjbQvLunjg6UXOetn4rv1xR8VisdzO/q7vcyTsU/cMuSXV6c+GLnrdeh2dKSrxtNJK6aHrhllxcDDg00s20eOyTh5aTadZm+WGXn0sLxDTO8jBPyytf5MPWovHgdgK+a2UaS3eHQQDNWXjzIy6xFpOr5OEz9/L1Pdwpwd74zK+jD1HPKi7PKjEHy9wDm5vtQdnnx1F5unUSl4xi/fipu76VnAZLrwnV+UpSNARd2FV1eXEwDzfHK+Uoevm40kbhy7fMleT2AkwB8BuBNABeY2QdtfaaU8uLrkjbQXKDcnsQMa5iBP+tISYk5xNuRkmVb+S4E8DMz20byWmROf2z7fHM10BSREPnI+Z5dc4pTzLmv7kFvOV+SpwI4w8zazKepgaakiiQ+ulHHPUuuXS+r9TKOa9ohz7njtdHzqmJciNzCtLzKGnxFRMrJNe2QvTGgNSQfB9Ajz49+bmYPRe/5OYBtAAo+yHDN+V4G4DvI5HdfRqa6bReoe7F44GuVIxLnc7eDmR3f1s9Jng9gNICh5vBrvsvZDr0B/D2AA8zsE5IzkCm2OADqXiwlIomJNWPTnoYE5rq66V7GKeNuhxHIxL9joiN4C3JNO+wAYGeSnyOz4l2HzNO8IdHP1b1YiubrRhOJK2N58S0AdgSwMIpzz5jZ99r6gMvBOmtJ3oBMt4pPACwwswUknboXZyeyp0yZggkT9HBFvkASv+mjIgvJ9cN6P3u/y7XVzMz6J/1Mwa1mJLsAmAXgbAAfALgfwEwAt5hZ56z3bTazLgXGUt5BRJz42Go2qt8op5jz6OpHgywvPh7A22a2CQBIPgDgG3DsXpxNOV+JI4ktl6nIQnJ1utFTkUXAMccl+K4GcERU5fYJMsdKLgPwERJ2L1bOV/LxdaOJxFV063gze5bkTAAvILN/7c/I7IfbDQm7F6uMVOIOWf0wLqo5Pe1pSGBuq5vlZZyQe7ipvFhEguQj5zu0zzCnmLOofkGQOV9vftVPT7Ul1xWrp+LXffXnQnL9eI2f3Q4hr3xVXiwiVaviO1m0Ul48CQm7FyvtICKufKQdvtl7qFPM+ePaReGlHdooLwYSdi8OeduHpIMktl59XtrTkMDsduXdXsaphrRDvvLivZNeTFvNJB9fN5pIXEUH3zbKi78Bh+7F2eXF1w8YgPN69fL6HyCVrceSJfjwkhPTnoYEZvf/eMTLOCH/tl1KefFCOHQvjo0V7v8JEQmKj5zv4F7HOMWc59YtCS/ni1bKi83snuY3tNW9OFvIfwtJOlReLPl4Ky+u5LQDWikvLqZ7sXK+ko/Ki6W9NFoZD5VMqJTy4luTdi9+oe+YkiYr1efQNXPwDzpMXWKu93TGc8i/bau8WESC5CPne1CPbzjFnBc3/FeQOV9vpvc4p5yXkwowdsM0/E6HqUvM9yvsMPViaOUrIkHysfI9sPsRTjHnlY3PhLnyJXkpMqXEBPB7M7uJ5B5I2L34mZ6nljRZqT5HrJ+NWd3HpT0NCczpG+/1Mk7IK9/tCr2B5IHIBN7BAA4CMJrkAGTOdlhkZgMALIpei4gEo9GanL7S4FJkcSaA4Wb2nej1PwL4FMBFAIZktRFabGYDC4wV7l9DIhIUH2mHv+16mFPM+Z9Ny4JMO7wC4BqSeyKzz3cUMuXE6l4sJSOJp7qrk4XkOmqjn04WIacdXI+UvAjADwFsBfAaMkH4AnUvFpH24mPlu+9ehzrFnDfffSHIlS/M7DYAtwEAyV8BqEcR3YsnajO9xFxXNx1zu+nPheQa3eCpyKIKVr7dzKyBZD8ACwB8HcAVAN4zs8kkJwHYw8wmFhgn3P8TIhIUHyvfmj3/zinm1L33UpgrXwCzopzv5wB+aGabSU5Gwu7FIZf6STpI4sPvDk97GhKY3afM9zJOyDFHRRYiEiQfK98+exzoFHPq338l2JWvF/O7nV3Oy0kFGN5wnw7WkRZ0sI7vi2nlKyKOfKx8e3Y+wCnmrP/gtTBXvq2UF1+FhN2L13/zmBKmKtWo5x+X4Op+OnBJcl25epqXcSp6t0NUXjwdmfLizwA8BuD7AM4FsDVJ92KtfEXElY+Vb9e/GegUczb978ogV777A3jGzD4GAJJLkOlckVjI+RdJB0lcVKMKN8l1W52nCreAY07Bg3WQKS8+muSeUSuhUQD6Rj+7mORLJG+PGm22QHICyWUkl9XW1nqatohIYU1mTl9pKKW8eDLUvVhE2omPtEOX3fo7xZzNW1eVPe2QeLdDc3mxmf0263t7A5hrZge29dmZPc9V8JUcZ26YhhUDRqY9DQnM/m/M8xJ8/2a3fZ1izv9ufdNL8CV5OYDrAXQ1s3fbfG8J5cU7NZ9qRvIyAIebWZsbNrXyFRFXPoLv7rt+xSnmfPjRWyVfi2RfALcC2A/AoELBt5Ty4ruTdi8OOfkt6SCJd0cfnfY0JDB7zV3qZZwyH5R+I4CJAB5yebOKLEQkSD5WvjvvXOMUc/7yl9XfRXTueKTWzJx3CJAcA2ComV1K8h0Ah/la+Xrx3hitcCTXnnOWYvPYIWlPQwLTZfpiL+O4Li6jQNtmsCX5OIAeeX70c2ROeRyWZG5a+YpIkHysfHfcqa9TzPn0L2uKvhbJryHTx/Lj6Ft9AKwDMNjMNrT2ubKufJXzlTiS+OS5mWlPQwKz8+AzvIxTjphjZi8D+GsbtSDTDmTZt9JJBfB1o4nEpVVA4aKsaQf5AskJSRL60jHoz0XH4VJeLO1DbZwlH/256CAUfEVEUqDgKyKSAgXf9CivJ/noz0UHoQduIiIp0MpXRCQFCr4iIilQ8C0zkiNIriS5iuSktOcjYYi6wTSQfCXtuUh5KPiWEcntAfwGwEgABwAYR/KAdGclgfgDgBFpT0LKR8G3vAYDWGVmb5nZZ8h0hT455TlJAMxsKYD3056HlI+Cb3n1BrAm63V99D0R6WAUfMsr38lC2usn0gEp+JZXPYC+Wa+bz/0UkQ5Gwbe8ngcwgOQ+JL8MYCyAOSnPSURSoOBbRma2DcDFAOYDWAFghpm9mu6sJAQk7wXwJwADSdaTvCjtOUn7UnmxiEgKtPIVEUmBgq+ISAoUfEVEUqDgKyKSAgVfEZEUKPiKiKRAwVdEJAX/B2vMC5kpc2rnAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Sample mar data -----------------------------------------\n", "X_miss_mar = produce_NA(X_complete_cont, p_miss=0.4, mecha=\"MAR\", p_obs=0.5)\n", "\n", "X_mar = X_miss_mar['X_incomp']\n", "R_mar = X_miss_mar['mask']\n", "\n", "print(\"Percentage of generated missing values: \", (R_mar.sum()).numpy()/np.prod(R_mar.size())*100, \" %\")\n", "\n", "ax = sns.heatmap(X_mar, mask=R_mar.numpy()==1, linewidths=0.005, linecolor='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### MNAR\n", "\n", "For generating MNAR data, an aditionnal argument should be given:\n", "\n", "* `opt`: it indicates how the missing-data mechanism is generated; using a logistic regression, logistic regression for generating a self-masked MNAR mechanism (\"selfmasked\") or quantile censorship (\"quantile\")\n", "\n", "#### MNAR with logistic model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With `opt = \"logistic\"`, missing not at random data are generated with a logistic masking model. It implements two mechanisms:\n", "* Missing probabilities are selected with a logistic model, taking all variables as inputs. Hence, values that are\n", " inputs can also be missing.\n", "* Variables are split into a set of intputs for a logistic model, and a set whose missing probabilities are\n", " determined by the logistic model. Then inputs are then masked MCAR (hence, missing values from the second set will depend on masked values).\n", " \n", "In either case, weights are random and the intercept is selected to attain the desired proportion of missing values.\n", "\n", "An additional argument is required: \n", "\n", "* `p_obs`: proportion of variables that will be used for the logistic masking model. Note that at the end, these variables are missing, since they are masked MCAR" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percentage of generated missing values: 39.5 %\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAAD7CAYAAADNT5fNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3debRU5Znv8e+jxqioF1GZBFEjC/WaxqgXk2a1ojgHhzii0cQhwU63aWM0iuauxJvENDEmaK8kRpxtB0RAgxgERYGOtx3AOSLXEThMxwHjGBHOc/+ofbRO1T5Vu6reqr2rzu+z1llSdfZ+9+tZVW899e73eR9zd0REpLE2SrsDIiI9kQZfEZEUaPAVEUmBBl8RkRRo8BURSYEGXxGRFNQ0+JrZ4Wa2xMxeMbPxoTolItLqrNp1vma2MfD/gEOANuBJ4BR3fzFc90REWtMmNZw7AnjF3V8DMLPJwDFAt4OvmSmjQ0QScXertY1P33ot0Zjzhe12qflalapl8N0BWJ73uA3Yr/AgMxsHjAOYOGIYZ+w6sIZLSqvZ5o5H+GjmxLS7IRmzxZjz0+5C3dUy+MZ9UhR9yrj7JGAS5CLf859YUsMlpRX1hDeapKRjQ9o96FYtg28bMDjv8SBgZakTtI+EFDIzvS6kiFmgWYAN68O0Uwe1DL5PAkPNbGdgBTAWOLXUCcH+oNJS9LqQenHvSLsL3ap68HX39WZ2LjAb2Bi40d3/WuqcV/7nodVeTlrUrn+dw8ePT0m7G5Ixm+93UpiGOrI7+Fa91Kyqi2m1g4gkFGK1w7rlzyYaczYdPLypVjtUTHN7UkhzvhIn2FRUi95wExHJtlac8wUwszeA94ENwHp337fM8bVcTlqUXhdSL96iqx06HejubyU5cMluhwW4nLSSYS/N5tH+x6XdDcmYkaunh2kowzfcGjrtMOyl2Y28nDSJYG80kUKtOu1ALqNtTrSK4doom62L/PTi3+49lDN2UXqxfK7P1Pk8veNRaXdDMuYry+4L01CGb7jVtNTMzAa6+0oz6ws8CHzf3ReUOF63tUUkkRBLzT5Z/EiiMeeLux/YXEvN3H1l9N92M7uH3E5n3Q6+S/c5qJbLSQsasuhh3j1jdNrdkIzpffPcMA214g03M+sFbOTu70f/PhT4Walzhix6uNrLSQsL9kYTKdSiN9z6AfdEy4Q2Ae5w9wdKnaDF9FLIzFj1Twek3Q3JmAH/NT9IO+7ZnfOtZW+H14DhlZyj9ZwSJ9QbTaRIC692qMjKkaMaeTlpAgMfncdHt/3vtLshGbPFab8I01CGpx20sY6IZFKI1Q5/X3RvojFns32Ozd5qBzO7ERgDtLv7ntFzfYC7gJ2AN4CT3H1tubYe6XtCLX2VFnRg+1TdC5Ai4TZT/zRMO3WQpHT8zcDhBc+NB+a6+1BgbvRYRCRbOjqS/aSgbOTr7gvMbKeCp48BRkX/vgWYB1xcrq0D26dW1DnpGXQjVuqmBW+49XP3VQDuvirKcIuVn178na1HcPAWu1Z5SWlFY1ffwcZfUMq5dLXh05LlIJPL8A23uq92KKxefP17T9T7ktJkgr3RRAq14OC7xswGRFHvAKA9yUlPDT66ystJq9p7+Qzev0CvC+lqq9/MCNKOZ/iGW7WD7wzg28CE6L9/SnLS3svD/EGltYR6o4kUaeY5XzO7k9zNte3MrA34KblBd4qZnQ0sA05McjEtKZJCZsa6N19NuxuSMZtu/6UwDTXztIO7n9LNryreikp3tSVOsDeaSKFmjnxDWjRIc3vS1T5tM3hx1yPS7oZkzB6vzArTUDNHviHt06a5PSkW7I0mUqiZI99u0osvA74LvBkddqm7/7lcW5rzlUJmxod/ODftbkjG9PqX34VpaH1zb6Z+M/A74NaC5ye6+5WVXExzvhIn2BtNpFAzR77dpBdXRZGvFDIzXhpWuHWI9HS7LSlZlyG5DM/5JtlYpzvnmtlzZnajmW3T3UFmNs7MFprZwkmTioobi4jUj3ck+0lBov18o8h3Zt6cbz/gLXKl438ODHD3sxK0o9BXRBIJsZ/vx/dMSDTmbP6N8dnbzzeOu6/p/LeZXQfMTHLe+CHdLRmWnmrC0ju5ceA30+6GZMxZK28P01Azz/nG6dzXIXr4DeCFJOdNWHpnNZeTFhfsjSZSKPBqBzPbGFgIrHD3MbW0VW168Sgz24vctMMbwDlJLqYbblLIzLi/78lpd0My5uvtd4VpKPyYcx6wGNi61oaqTS++oZqLaamZxAn2RhMpFHC1g5kNAr4OXA78sNb2Gprh9t73tKRIutr6mgf0jUiKBAvUEg6++UUfIpOivcjzXQVcBGwVomsNHXy3vibQ2j1pKfpGJHWT8IZbftGHOGbWmeW7yMxGhehakjnfweSy2/oDHeQ+Ea6upoKxqhdLoQPbpzJLc75S4IhQU1EbNoRpB0YCR5vZkcBmwNZmdpu7n1Ztg2XX+UaVKga4+1NmthWwCDgWOAN4x90nmNl4YBt3L1lEU+t8RSSpIOt8b7oo2TrfM69IfK0o8r2w7qsdoiVlncUy3zezxcAOVFHB+N931HpO6eqSZbfzK70upMDFywItP2yV9OIo0+0rwOMUVDAGYisY56cXP/HBy7X1VkSkEnVIL3b3ebVGvZAwvRjAzLYE5gOXu/t0M3vX3Xvn/X6tu3e7x0N0jKYdRCSRENMOH006P9GYs8W4idlMLzazLwDTgNvdfXr0dMUVjLWkSAqZmV4XUqTRS83SkGS1g5FLqljs7r/N+1XFFYy1pEji6HUhdRNutUNwSSLfkcDpwPNm9kz03KVUUcFYEY4UMjOm9dOGS9LV8WsC7QPTzJGvu/8F6C40qaiCsSIciRPsjSZSqJkH35AU+UohM+PMIcen3Q3JmJuWTgvTUIbHnIYOvop8JU6wN5pIoWaOfEukF19GhRWMFflKITPj6MFfT7sbkjEzlt8fpqGO7I45SSLf9cAF+enFZvZg9LuKKhgr8pU4wd5oIoWaebVDifTiiinylUJa5ytxQgVq3szTDvkK0otHkqtg/C1yZTUuiNvVLGafTJEu9I1I6qbJpx2Az9KLpwE/cPf3zOwacpWLOysY/wYoqmBcsE9mdv8SkgpFvhIn2AdysxfQjEsvrraCsYhIwzRz5NtdenE1FYz19VLi6HUhdbO+iW+40X168SmVVjB+efdDq+ymtKqhi+fw5A7HpN0NyZj/taLsVjHJNPO0Q4n04pJreuMMXTyn0lOkBwj2RhMp1MzTDiEt2e2wRl5OmsCwl2brhpsU0VIzwMw2AxYAX4yOn+ruPzWznYHJQB/gKeB0d19Xqq1hL82uvcfScjTnK3XT5JHvJ8BB7v5BtOrhL2Y2C/ghuQy3yWb2R+Bs4JpSDSnCkUJmxi0DVcNNuvr2ylA13LI75iSZ83Xgg+jhF6IfBw4CTo2evwW4jDKDryIciRPsjSZSqJnTiwHMbGNyJeN3BX4PvAq86+7ro0PaSJByrMhXCpkZNyrylQJnBfpA9maOfAHcfQOwl5n1Bu4Bdo87LO5cpRdLOaHeaCJFmn3w7eTu75rZPOCrQG8z2ySKfgcBK7s5R+nF0i3N+UqccHO+zb3aYXvg02jg3Rw4GPgV8AhwArkVDyqgKVXTnK/UTZNHvgOAW6J5342AKe4+08xeBCab2S+Ap8mlIJf0xMBja+qstJ4RK+/l9eEHp90NyZidn30oTEPNPPi6+3PktpEsfP41YEQ9OiUiEoJvaOJph5BGrLy3kZeTJhEsyhEp1MyRb0haaiaFtJ+vxAmXXpzd11Yt6cU3AwcAf4sOPcPdn4lv5bO2auuttCS9LqRumnnwpfv0YoAfufvUpBeb2/fEavooLWx0+938ckctNZOuLl0WaqlZmGbqoZb04oqNbr+7mtOkxQV7o4kU8PXZHX2rSi9298fN7HvA5Wb2E2AuMN7dPynVztJ9Dqq1v9Jihix6WHO+UiTYVFR2x97q0ovNbE/gEmA1sCm5DLaLgZ8Vnqv0YilHc75SL019wy1fXnrx4e5+ZfT0J2Z2E3BhN+d8ll782MDjsvuXkFR8bdU9inyliCJfuk8v7iygGRXYPJYEBTS/tuqemjssrUeRr9RLqMjXzAYDtwL9yQ3pk9z96lrarCW9+OFoYDbgGeCfyzX0aP/jaumrtKCRq6cr8pUiGYx81wMXuPtTZrYVsMjMHnT3F6ttsJb04orvno1cPb3SU6QHUOQr9fLZjuO1tuO+ClgV/ft9M1tMbg/z+g2+IU3pf2r5g6RHOWn1HVw1WOt8pasfLA+1mXqQZrows53IBaSP19LORiE6IyKSSR3JfsxsnJktzPuJXaFlZlsC04AfuPt7tXTNks63RXO+C4EV7j6mmurFZqbJPRFJxN1rno9685ADEo052z84v+y1ogzfmcBsd/9trX2rZNrhPGAxsHX0+FeoerHUSBvrSJxwG+sEaYZoVdcNwOIQAy8kz3AbBHwduBz4YdQRVS+WIPS6kHrxDcFeWyOB04HnzaxzA7FL3f3P1TaYNPK9CrgI2Cp6vC0JqxfnZ7hde+21jBunZDf5nCJfiZO1yNfd/0JuWW0wSZIsxgDt7r7IzEZ1Ph1zaOw7KD/Dzcz8nHPOqbKr0qoU+Uq9eEd2X1tJIt+RwNFmdiSwGbk536tIWL04391aaiYFTlx9B2tPGpV2NyRjtpkyL0g79VhqFkri1Q4AUeR7YbTa4W5gWt4Nt+fc/Q9lztf3SxFJJMRqhxVfOyjRmLPDfz/c8BC5liSLi6mwevGHV5xZw+WkFfW66CY+nKj7ANJVr/MnBWmnZSLfmi+myFdEEgoR+S7bd3SiMWfHhXObKvKt2LJ9tZm6dLXjwof1jUiK9LropiDtNPsNNxGRppTlwbeW9OKbqbx6saYdRCSRENMOrw8/JNGYs/OzD2Z62qEwvRgqrF6sxfRSyMz4tyEnp90NyZj/WHpXkHayHPlWlV5c7cW0mF7ihHqjiRQKEDzXTbXpxZ3KVi9WerGUYmbcNkD7+UpXp60Ks5/vhnB7OwRXds43Si8+0t3/pSDJYgBdqxe/6u5F1YsL2tK8g4gkEmLOd8luRyQac4a9NCuTc75F6cVmdpu7nxb9vmT14nxrDtq/+p5KS+r38ALdC5Ai4TbWyW7km6SG2yXAJdAlvfi0aqoX93t4QY3dlVakewFSL1n+XK9lne/tlVYvXvVPB9RwOWlFA/5rviJfKdITIl+lF4tIJoWY831+56MSjTlffv2+TM75BvPut5ReLF31vvVhFvQ7Pu1uSMbsv2ZakHay/KVKka+IZFKIyPeZIUcnGnP2Wjojm5Gvmb0BvA9sANa7+75m1ge4C9gJeAM4yd3Xlmrn1T0PqaWv0oK+9MKDvLjrEWl3QzJmj1dmBWkny0kWG1Vw7IHuvpe77xs9Hg/MdfehREkWwXsnIlID92Q/aUg07RBFvvu6+1t5zy0BRkXLzQYA89x9WJl2NO0gIomEmHZYOOjYRGPOvm33ZnPagVxxzDnR4HltVBSzn7uvAogG4L5xJ+anF/9k2y9z4tY7Bui2tIovv34/Hy/6U9rdkIzZfJ9jgrSzoaOSL/eNlTTyHejuK6MB9kHg+8AMd++dd8xad9+mTDuKfEUkkRCR72MDj0s05nx15fRsRr7uvjL6b7uZ3QOMANbkZbkNANoTtFNTZ6X1mBnfGnJc2t2QjLl16fQg7XRk+IZb2cHXzHoBG7n7+9G/DwV+BswAvg1MiP5b9ruj0kglTqg3mkihLK92SBL59gPuiQbOTYA73P0BM3sSmGJmZwPLgBPLNaTIVwqZmV4XUiRUoJbh4sWJNtZ5DRge8/zbwOhKLqbIV+LodSH14mT3tdXQ9OJ7+53SyMtJEzh2zZ3M7qsyQtLVYe1hqpusz/C0g9KLRSSTQqx2mNvv5ERjzug1d2VztUM36cWXAd8F3owOu9Td/1yqnQ8nnFF1R6U19Rp/s+Z8pYjmfLs6MD/DLTLR3a9M2kCv8TdXcDnpKTTnK/WiOd/I++cf1cjLSRPYauJ9inylSE+IfJPm3nWmFy+K0oU7nWtmz5nZjWYWm91mZuPMbKGZLbzxuaU1d1hEJKkNWKKfNNSSXrwEeIvcwPxzYIC7n1WmHYU4IpJIiBtu9/U/JdGYc9TqO7N5wy0uvdjdP6uGaWbXATPLtfPm4apeLF1t/8ACLhoyNu1uSMZcsXRykHY6Aka1ZnY4cDWwMXC9u0+opb2q04s793WIDvsGCaoXb/+AqhdLsVBvNJFCob5qm9nGwO+BQ4A24Ekzm+HuL1bbZi3pxf9pZnuR+/97AzinXEOvDz+42n5Ki9r52Yd0w02KZPCG2wjglSjjFzObDBwD1G/wLZFefHqlF9v52YcqPUV6AC01k3rpCPfa2gFYnve4DdivlgYbutTstX9Q5Ctd7fLcQ7x/4bFpd0MyZqsr7w3SzoaEx+UXfYhMiopGfHZIzGk1fWVr6OC7y3OKfKVYqDeaSKGOhIFvNNBOKnFIGzA47/EgYGXVHSN5enFv4HpgT3Kj/VnklppVVL34rTFa7SBdbTdzgeZ8pUi4Od9g0w5PAkPNbGdgBTAWOLWWBpNGvlcDD7j7CWa2KbAFcCm56sUTzGw8uerFF5dqZLuZWu0gxTTnK/US6mPd3deb2bnAbHJLzW5097/W0maSpWZbA/sDZ0SdWAesM7NjgFHRYbcA8ygz+GrrQCl0WPtdinylSLDIN+DnerRxWMnNwyqRJPLdhdzOZTeZ2XBgEXAeVVQvFomjyFfqJct7OyQZfDcB9ga+7+6Pm9nV5KYYEsmfyF574iiFONJFn6nzFflKkVAfyBsy/LmeZGOdNqDN3R+PHk8lNxiviaoWk7R6sYhII3Uk/ElDkiSL1Wa23MyGufsScnXbXox+Kqpe3Gfq/Bq7K61I0w5SL80+7QC5Xcxuj1Y6vAacSS5qVvViqYmZMWmHb6bdDcmYcStuD9JOhku4Jd7V7Blg35hfqXqx1CzUG02kUCtEvkG8c9wBjbycNIE+0+dzvSJfKfCdQB/ISdOL09DQwbfPdM35SrFQbzSRQiHX+YZWS3rxYVRYvXjlyFFVd1Ra08BH5/HOCfpGJF2FujnfCtMOcenFh1Fh9eKBj86rvIfS8rQKRuqlqQffEunFFV/s6R1VvVi6+sqy+7hysOZ8pasLlwda7RCklfqoJb0YctWLvwUsBC6I29VM6cVSTqg3mkihLM/5lq1ebGb7Ao8BI/PSi98DfkeF1YvJ9geRpEDrfCXOuBW3B6le/O9DTks05lyy9LaGD9NVpxe7+xp33+DuHcB15GociYhkRgee6CcNVacXV1O9WEkWEkdJFlIvTX3DLRKXXvwflVYvVnqxFDIzrhmkaQfp6nttuuEGdJteXHH1YkW+EifUG02kUCtEvkEsHnpEIy8nTWD3l2cp7VyKhMqGXW/ZjX2TrPMdRq5QZqddgJ8At1JhAc3dX55VbT+lhSntXOolu0NvshtuS4C9AMxsY3KVO+8hV82iogKamvOVQmam14UUCVe9OLsqnXYYDbzq7kurKaCpOV+Jo9eF1Etay8iSqHTwHQvcGf07UQHNfIpwpJAiX4kT6gM5y6+sxINvtMzsaOCSSi6g9GIpR5Gv1EurTDscATzl7muix2s6Ey1KFdDMr17ctt9BWf4gkhQMfuIRrYKRIqFuzm/IcOxbyeB7Cp9POQDMoMICmoOfeKSizknPoFUwUi9NH/ma2RbAIXTNYptAhQU0P/g/p1bTR2lhW/70Ds35SpFwc77ZfW0lzXD7CNi24Lm3qbCApohIIzV95BvKlj+9o5GXkyahG25SL6201Kwm+nophcyMS4ecknY3JGN+ufTO8gclkOURp5b04t5UWEBTEY7ECfVGEym0PsPDby3pxWdSYQFNRb5SSEkWEkc33IrlpxdXfDFFvhJHrwupl1a64ZafXgwJCmjmWzt2VIWXk1a3zeR5inyliCLfPDHpxdeQK5zZWUDzN0BRAU2lF0s5inylXlol8u2SXpyXZoyZXQfMjDspP72YbN98lBRozlfihPpA3pDh11bV6cUqoCmh6HUh9dKodb5m9mvgKGAd8Cpwpru/W+qcWtKLr6i0gOacvicluZz0IIe2T1HkK0WacM73QeASd19vZr8iNz1bcn/zWtKLKy6geWj7lEpPkR5Aka/US6PmfN19Tt7Dx4ATyp2jDDdJlZnx4UTdj5Wuep0/qfxBCSSddohZGDApul9VjbPompgWq6GDr4hIIyWddihYGBDLzB4C+sf86sfu/qfomB8D64Hby10z6Zzv+cB3yM3vPk8uu20AMBnoAzwFnO7u68q0k+Ry0sOEinJECoVc7eDuB5f6vZl9GxgDjPYEX/OT7O2wA/BvwB7u/rGZTSGXbHEkufTiyWb2R+Bscmt/S3W+3OWkhzEzLhoyNu1uSMZcsXRykHYauNrhcHI32A6I7pGVlXTaYRNgczP7FNgCWAUcBHTujn4LcBllBl9FvhIn1BtNpFADkyx+B3wReDAa5x5z938udUKSjXVWmNmV5KpVfAzMARYB77r7+uiwNmCHuPPzJ7KvvfZaxo3TzRX5nJIsJE6zLTVz910rPSfJtMM2wDHAzsC7wN3kst2Krt9Npz6byDYzP+ecssuBpYfRNyKpl2bfTP1g4HV3fxPAzKYD/wj0NrNNouh3ELCyXEPLRxxYS1+lBQ1+4hHeP/+otLshGbPVxPuCtJPlb1VJBt9lwFejLLePyW0ruRB4hNxC4smoerHUINQbTaRQU5eOd/fHzWwqueVk64GnyU0j3A9MNrNfRM/dUK6tp3dUhCNdfWXZfZw95Pi0uyEZc8PSaUHayfK0gzUyLDez7P4lRCRT3L3mmwGjBx2aaMyZ2zan4TcelF4sqdJqB4kT6iZsliNfpReLSMtq+koW3aQX/xE4APhbdNgZ7v5MmXaq76m0LL0upF6aejP1EunFAD9y96lJL6avl1JI0w4SR9MOXY/LTy8uu6Y3jiIciaPXhdRLUw++cenF7j7HzE4FLjeznwBzgfHu/knh+fnpxb8eOpTTBw4M+j8gza3//PmKfKVIsPTiDL+2yi41i9KLpwEn83l68VRyA+5qYFNy635fdfeflWkru38JEcmUEEvNRgw8INGY88TK+ZlcahabXuzut0W//8TMbgIuLNfQs0PGVN1RaU3Dl87MdHQi6Wi2jXWqUXV6cWf1Ysv9lY4lQfXi4Utjq8tLD6c5X6mXDd7ATSUrVEt68Swz2x4w4Bmg5N6VAE8NPrq23krL2Xv5DH6kzdSlwK8D7fGc5W9VSi8WkUwKMec7vP8/Jhpznl39fzM55xvM5P6nlj9IepSxq+/gmkHfTLsbkjHfaytbfzKRLM/5KvIVkUwKEfnu2e+ricacF9Y8ls3I18zOA75Lbn73One/ysz6kKtNvxPwBnCSu68t1U6W518kHWbGtH6npN0NyZjj19wZpJ0sR74blTvAzPYkN/COAIYDY8xsKDAemOvuQ4mSLOrZURGRSm3wjkQ/aUgS+e5OrhLnRwBmNh/4Brm6bqOiY24B5pErndwtLSmSOKGiHJFCHRn+tp1k8H2BXBrxtuTW+R5JroxQP3dfBRCt9+0bd7KqF0sp2lhH4vSEJItEN9zM7GzgX4EPgBfJDcJnunvvvGPWuvs2ZdrJ7l9CRDIlxA23L223d6Ix59W3nsrmDTd3v4GoRpuZ/RJoA9bkZbkNANrLtXORFtNLgSuWTmZmX70upKsx7YGSLFog8u3r7u1mtiMwB/gacCnwtrtPMLPxQB93v6hMO9n9S4hIpoSIfIds+w+Jxpylbz+XzcgXmBbN+X4K/Ku7rzWzCcCUaEpiGXBi2Ua0pEgKHL/mTs35ShFtKRn6Yop8RSShEJHvoD57Jhpz2t55IbORbxCz+57cyMtJEzis/S5trCNFtLFO6Isp8hWRhEJEvgN675FozFn17ovZjHy7SS++LHruzeiwS939z6XayfKnkKRD63wlTk9Y55ukenF+evE64AEzuz/69UR3vzLpxZThJnH0upB6aerN1Ok+vbhiinCkkCJfidMTVjuU3ViHXHrx/ma2bVRK6EhgcPS7c83sOTO7MSq0WcTMxpnZQjNbOGnSpEDdFhEpr8M90U8aakkvngC8BTjwc2CAu59Vpp3sfgyJSKaEuOG2zZa7Jhpz1n7wSsPnvipe7dCZXuzuf8h7bidgprvvWercqQO+qcFXujhx9R0sHnpE2t2QjNn95VlBBt//seWXEo05f/vg1SCDr5ldCPwa2N7d3yp5bA3pxZt17mpmZucD+7l7yQWbinxFJKkQg+/WvXZJNOa89+FrNV/LzAYD1wO7AfuUG3xrSS/+TzPbi9y0wxvAOeUaUYQjhXZ/eRZvjdk/7W5Ixmw3c0GQdhq82mEicBHwpyQHK8lCRDIpROS7+eZDEo05f//7snOI9h2PTHL3xCsEzOxoYLS7n2dmbwD7hop8g3j7aEU40tW2MxawduyotLshGbPN5HlB2kkaXEYDbcnB1sweAvrH/OrH5HZ5PLSSvinyFZFMChH5fnGzwYnGnE/+vrzqa5nZl8nVsfwoemoQsBIY4e6ruzuvoZHvhxPOaOTlpAn0Gn8zHz8xNe1uSMZsPuKEIO00Irh09+eBz8qoJZ12UOQrIpkUIvLdZNMdEo0569etCLbON5ODr3zOzMZVMqEvPYNeFz1HkvRiqQ+VcZY4el30EBp8RURSoMFXRCQFGnzTo3k9iaPXRQ+hG24iIilQ5CsikgINviIiKdDg22BmdriZLTGzV8xsfNr9kWyIqsG0m9kLafdFGkODbwOZ2cbA74EjgD2AU8xsj3R7JRlxM3B42p2QxtHg21gjgFfc/TV3XwdMBo5JuU+SAe6+AHgn7X5I42jwbawdgOV5j9ui50Skh9Hg21hxm3dorZ9ID6TBt7HagMF5jzv3/RSRHkaDb2M9CQw1s53NbFNgLDAj5T6JSAo0+DaQu68HzgVmA4uBKe7+13R7JVlgZncC/w0MM7M2Mzs77T5JfSm9WEQkBYp8RURSoMFXRCQFGnxFRFKgwVdEJAUafEVEUqDBV0QkBRp8Rfm1k/4AAAAJSURBVERS8P8BYQRL7lg7zk0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Sample mnar data with logistic model -----------------------------------------\n", "X_miss_mnar = produce_NA(X_complete_cont, p_miss=0.4, mecha=\"MNAR\", opt=\"logistic\", p_obs=0.5)\n", "\n", "X_mnar = X_miss_mnar['X_incomp']\n", "R_mnar = X_miss_mnar['mask']\n", "\n", "print(\"Percentage of generated missing values: \", (R_mnar.sum()).numpy()/np.prod(R_mnar.size())*100, \" %\")\n", "\n", "ax = sns.heatmap(X_mnar, mask=R_mnar.numpy()==1, linewidths=0.005, linecolor='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Self-masked MNAR" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For self-masked MNAR values, the missingness of the variable $X_j$ only depends on the values of $X_j$. \n", "\n", "With `opt = \"selfmasked\"`, self-masked missing not at random are generated with a logistic self-masking model. \n", "\n", "Variables have missing values probabilities\n", " given by a logistic model, taking the same variable as input (hence, missingness is independent from one variable\n", " to another). The intercepts are selected to attain the desired missing rate." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percentage of generated missing values: 39.0 %\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAD7CAYAAABHYA6MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZhcVbnv8e8LEQFlSiAhEwkI4gCSI7k5cPCSMMkUgUCABHAAFPUKh+OABnxEDhwUh6NyRJEZlCYBAolhHoKQq1eQgAgCxgMhQ6ehAwEhXMYk7/lj74bq6t1du2qvqtq16/d5nn7SXb1r7ZV+dq16a+31rtfcHRERyZ/1mt0BERFJpgFaRCSnNECLiOSUBmgRkZzSAC0iklMaoEVEcirTAG1mB5jZIjN7ysxmhOqUiIiA1boO2szWB/4O7Ad0Ag8C0939iXDdExFpX4MyPHcC8JS7LwYws1nAoUC/A7SZKStGRFJxd8vaxtsvLE415rxny+0yn6sesgzQI4HlJT93Av9cfpCZnQScBPD9MTty7FYjM5xSimabhffw2jVnNrsbkjMbH3N2s7uQC1kG6KR3nD7vVu5+MXAxRBH06UsXZTilFJFejFI369Y2uweZZBmgO4HRJT+PAroGesIrp07OcDopok3Pv5nlE/ZqdjckZ0b/6XdhGlq7Jkw7TZLlJuEgopuE+wAriG4SHuPujw/wHM1Bi0gqIeag3+p6PNWYs8GIjxZrDtrd15jZycAdwPrA5QMNzgC3DT261tNJQR248lpeu++KZndDcmbjiceHaWjdujDtNEnNEXRNJ1MELSIpBYmgl/8lXQQ9epdiRdC1uHvoUY08nbSAfVdex4tTJza7G5Izg2ffF6ahNr5JKCKSb97GUxxmtgRYDawF1rj7+ArHa4pDRFIJMcXx5tP3pxpz3vuB3Qo7xbGXu7+Q5sB5w6YFOJ0UySHds/j5qGOb3Q3JmVM6O8I01M43CeMIenzaAVoRtIikFSSC/vvv00XQH/xEISNoB+6MB96L4qzBXkpTvc8d+SGmDxmV8ZRSJNs9ejdXjVAELb19titUBN3aNwmzRtAj3L3LzIYCdwGnuPuCAY5XBC0iqQSJoJ/8XboI+sN7FS+Cdveu+N+VZjaHaIe7fgfo+4ZNzXI6KaCJ3bNZsfukZndDcmbkH+8N01CLp3rXPECb2fuA9dx9dfz9J4EBd72Z2D271tNJgQV7MYqUa/GbhFki6GHAHDPraecad799oCcogpZyE7tn08hsVmkN8biSmXtrz0Fn2YtjMbBLNc9RBC1JQr0YRfpo8USVhmYS/mnEYY08nbSACV1zFUFLH8HetNt4iqNqE7rmNvJ00iIUQUvdFD2CNrPLgcnASnffKX5sMHAtMBZYAhzl7i9VautHo7XeVXo7bXkHD48+pNndkJz5+PJ5YRpa+3aYdppkvRTHXAkcUPbYDGC+u+8AzI9/FhHJl3Xr0n3lVKpEFTMbC9xcEkEvAia5+7NmNhy41913TNGOJhtFJJUQiSpv/HFmqjFnw92n53KerdY56GHu/ixAPEgP7e/A0lTvvQePZ6dNPlDjKaWI/mvptbpJKH3oJmGk7jcJy6t63/PiwnqfUlqMbhJK3bTpAN1tZsNLpjhWpnmSIiUpZ2a8eLgqqkhvg28MU1HFW/wmYa0D9Dzgs8B58b+/TfMkRUqSJNSLUaSPNlhmNxOYBGxpZp3Ad4kG5uvM7ERgGXBkmpMpgpZyZsZbzz/d7G5IzmywVaB7VUWf4nD36f38ap9qT6YIWpIEezGKlCt6BB3S5dqYXcqc0NWhT1bSh1ZxRBo6QJ8QqkqCFIo+WUndFD2C7ifV+yzgC8Dz8WFnuPutldpSpCTlzEzXhfQR7E17TfE37L8SuAD4ddnjP3X3H1dzMkVKkkTXhdRN0SNod18Qp3pndsfQo0M0IwWy/8prmTusv/vQ0q4O654ZpqHAc9Bmtj6wEFjh7pPLfvc54EfAivihC9z90iznS7NZUn9ONrNHzexyM9uiv4PM7CQzW2hmC299XcupRKSBfF26r/ROBZ4c4PfXuvu4+CvT4Ay13yS8EDgH8Pjf/wROSDqwPNX7/NVK9ZbegkVLIuUCRtBmNgo4GDgX+FqwhgdQ0wDt7t0935vZJcDNKZ9Xy+mkwHSTUJIEuy+RMjou3dQtdnEcXJb6GfBNYJMBmjrCzPYE/g581d2XV9HbPmoaoHv24Yh/nAL8NeXzajmdFJyuC6mblKs4Sj/pJzGznpVsD5nZpH4OuwmY6e5vmtmXgKuAvavrcG+1pnpPMrNxRFMcS4AvpjnZFUpUkTLHd3XwnTHHNLsbkjPnLL0mTEPhPp3tARxiZgcBGwKbmtnV7n7cu6fyVSXHXwL8IOtJU23YH4o27BeRtEJs2P/6zO+mGnM2mv7vqc8VR9DfSFjF8c7MgplNAb7l7rtV0d0+GppJ+Pwn92zk6aQFbHXnAs1BSx+tkuptZmcDC919HvCvZnYIsAZ4Efhc5vYVQYtIHgWJoK/+droI+rhzc3kjJM0c9GiiLMKtgXVEdzfPr6Wyt6p6S7nTlmuzJOkrWAS9dm2YdpokzRTHGuDr7v6wmW0CPGRmdxGF7/Pd/Twzm0FU2ftbAzV02nJtliR9aRWH1E3Rd7OLJ717CsSuNrMngZHAoUSrOyBaTnIvFQboA0cfmKGrUkS3Lb+Ng0cf1OxuSM7csrzi3mvptPgAXVWqd7wnxz8BD1BW2RtIrOxdmuq97NVl2XorIlKN8KneDZV6FYeZvR+4Afg3d38l7cfS8lTvx196vJZ+SoEFi5ZEyvi61r6/kWqANrP3EA3OHe5+Y/xw1ZW97xs2tfaeSiFN7J6tm4TSR6sss6u3NKs4DLgMeNLdf1Lyq6ore0/snl1jN6XIdJNQ6qYNVnHsAXwaeMzMHokfO4MaKnsrUpJy2ixJkiiCjqRZxfF7oL+/VlWVvRUpSRJdF1I3RR+gQ1JVbyl3QlcH2wzeudndkJxZ9uJjYRpq8U9nquotTRfsxShSrugR9ACp3mdRZWVvzTVKOc1BS5Jwc9CtfW1lSfWGKit7a65Rkui6kLop+iqOAVK9q/arkZqDlt6+tKKD/Ubt3+xuSM7c1XlHkHa86FMcpcpSvfcgquz9GaIy5F9P2s0uodaXSC+hXowifbTBFAeQmOqdqrJ3Wa2v1v5rSXCag5YkjS4am1c1p3rXWtlbRKRhih5B95fqXUtlb90MkiS6LqRu1hT8JiH9p3pPr7ay9y1Dj66xm1JUB6+8lkt181jKfH5FoJyJFp/iUE1CEcmlEDUJ//+3j0w15rzv3Otz+TGuoZmENw+d1sjTSQuYvHIWs7Y+ptndkJyZ9tw1Qdop/DI7M9sQWAC8Nz5+trt/18y2BWYBg4GHgU+7+1sDtTV55azsPZbCCfViFOmj6DcJgTeBvd391Xg1x+/N7Dbga0SZhLPM7FfAicCFAzWk5VRSTsvsJIlSvSNpMgkdeDX+8T3xlwN7Az2fTa8CzqLCAK279ZJE14XUTdFTvQHMbH3gIWB74BfA08A/3H1NfEgnKdK/FSlJOUXQkiTUm3Zb1CR097XAODPbHJgDfDjpsKTnKtVbKlEELXXTDgN0D3f/h5ndC+wGbG5mg+IoehTQ1c9zlOot/TIzThqjYsLS28VLA9UvbYNVHFsBb8eD80bAvsAPgN8BU4lWcqQqGqtISZIEezGKlGuDCHo4cFU8D70ecJ2732xmTwCzzOw/gD8TpYMPSNuNSrkvrejg7qFHNbsbkjP7rrwuTENFH6Dd/VGiLUbLH18MTKhHp0REQvC1rT3FoVRvEcmlEKner5y4X6oxZ9PL7srl/GtDU721nErKmRnrv2dEs7shObP27cQ1B1ULtcyuv4zqsmPeS1S/dVdgFXC0uy/Jct4sqd5XAhOBl+NDP+fujyS38k5bWfoqBRXqxSjSR7g56MSMane/v+SYE4GX3H17M5tGtJgi0xaeWVK9AU5z99S34BVBSzklqkiScKneYZoZIKO61KFEGdUAs4ELzMw8wwWeJdW7aoqgJYmuC6kXX5NuhE5IqLs4zuEoPaZXRrW7P1DWzEhgOYC7rzGzl4EhwAu19b7GVG93f8DMvgyca2ZnAvOBGe7+5kDt3DdMCQnS28Tu2YqgpY9GR9BlCXX9HdMro9rMdnL30kpSSZ3OdHHXlOptZjsBpwPPARsQ/ce+BZxd/lyleksliqClXuqxF0dJRvUB9C711wmMBjrNbBCwGfBilnPVmup9gLv/OH74TTO7AvhGP895553pl6OPU6gkvXyls4NXz9AnK+nt/d8LleodppkBMqpLzSPKqv4jUZb1PVnmnyFDqndP0di4qOxhpCga+5XOQHXGpFCCvRhFygSMoPvLqD4bWOju84iyqX9jZk8RRc6ZS0hlSfW+Jx68DXgE+FKlhn4+Sqne0tspnR08tu3Bze6G5MzOz9wSpqFwqzj6y6g+s+T7N4Ajw5wxokxCEcmlEJmEqw6emGrMGXLLfbm8EaJMQmkqM2PKNp9qdjckZ+YsuylIO97aW3E0doAWEWmodhmg4znohcAKd59cS1VvLaeSJKGiJZFy7RRBnwo8CWwa//wDqqzqffkI3SSU3k7o6uDsbY6pfKC0lTOXXROknVYfoFPdJDSzUUSVu88FvgZ8Cnge2DpOadwdOMvd96/QjiahRSSVEDcJuydNSjXmDLv33lx+vE8bQf8M+CawSfzzEFJW9S7NJDxxswnss/H2tfdWCueYZ6/RzWPpI1xV7yDNNE2aRJXJwEp3f8jMJvU8nHBo4qusNJPQzPyyl/9UY1elqHRvQurF17X2tZUmgt4DOMTMDgI2JJqD/hkpq3qXUqQk5bTdqCRRBB1Js93o6UQbIxFH0N9w92PN7HpU1VsC0HUh9RJgGrupsqyD/hZVVvVe/VUlJEhvm/z0JlZ/c0qzuyE5s8kP5wRpp9UjaKV6i0guhVjFsWz8PqnGnG0Wzs9lqN3QTMIFw45o5OmkBezZfYPmoKWPcHPQuRx3U1Oqt4gUVtsM0Amp3ldSZVXvPbtvqLWfUmC6SSj10uofzrKkekOVVb3nDJtexemkHUzpnsnOw3ZrdjckZx7rvj9IO20RQcep3gfzbqp3TaZ0z6z1qVJgoV6MIuXaZZldeap3j4pVvUtTvU/fbBxT3je29t5K4UzomqubhNJHqGmvtWsLPkD3k+oNKat6l6d6f//lAaeppQ1pDlrqpR0i6D6p3mZ2tbsfF/9+wKrepf4yZnLtPZVC2mXpzdygexNS5ohA06GtPgddVaJKSar35LKq3j8F3nD3GRWer8+yIpJKiESVJ3c4KNWY8+H/vjWXI3mWddAd1Vb1fnDkoRlOJ0X0v1b8lvuGTW12NyRnJnanXhw2oLaKoDOfTBG0iKQUIoJ+bNtPpRpzdn7mplyO5A3NJOzcba9Gnk5awKj7f8dPRqsUmvT2teUdQdpp9QVCiqBFJJdCRNCPjDkk1Zgzbum81o2gzWwJsBpYC6xx9/FmNhi4FhgLLAGOcveXBmrnjqFHZ+mrFND+K6/VKg7pI9gqjhZfZrdeFcfu5e7j3H18/PMMYL6770CcqBK8dyIiGbin+8qrtFW9lwDj3f2FkscWAZPipXbDgXvdfccK7eT4TyEieRJiimPhqMNSjTnjO+fmMtROe5PQgTvjAfaiODtwmLs/CxAP0kOTnlia6v35TSewr6p6S4lpz6mqt/QVLNV7XTWTBPmTdoDew9274kH4LjP7W9oTlKd6X/qKqnpLb0r1lnpp9bf+VAO0u3fF/640sznABKC7JJtwOLCyUjtXD9dyKuntuGc7GLnFR5vdDcmZFS89HqSddS1+kzDNZknvA9Zz99Xx958k2hRpHlE17/NIWdX7uGfDrG2UYgn1YhQpF3IVh5ldDvRsHrdTwu8nEY2Dz8QP3ejufTaQq0aaCHoYMCf+GDoIuMbdbzezB4HrzOxEYBlwZKWGHh59SJa+SgF9fPk8zUFLH6GmvQIX9b4SuAD49QDH/F93D7YrXMUB2t0XA7skPL4K2Keak318+bxqDpc2oTloqRcn3LXl7gvMbGywBlNoaKq3IiUpZ2a6LqSPUG/aa1JOcZSuNotdHC9wqNbuZvYXoIto589M83cNHaAVKUkSXRdSL2kj6NLVZhk8DIxx91fj/fPnAjtkaTBLqvdZwBeA5+PDznD3Wwdq55VTDq69p1JIm/78Fl679t+b3Q3JmY2P/m6QdgLPQQ/I3V8p+f5WM/ulmW1ZmuBXrWoi6L0STvRTd/9x2gY2/fktVZxO2kWoF6NIuZBz0JWY2dZAt7u7mU0g2kpjVZY2GzrFseqwPRt5OmkBQ+Yu4LXrz2l2NyRnNj7yO0HaCRlBm9lMYBKwpZl1At8F3gPg7r8CpgJfNrM1wOvANM94gyVLqjfAyWb2GWAh8PWk3exKJ99/Mm4HPjt2eJb+ioiktjbsKo4Bt1109wuIluEFk3azpBGlqd7AKcAi4AWiwfscYLi7n1ChHd2uF5FUQmyWdNPW01ONOZ96bmYu71TXnOrt7gt6fm9mlwA3V2rn8e0OqrWfUlAfXXyrltlJH+ESVXI57qZWc6p3zz4c8WFTgL9Wauujiwdc5CFtSsvspF5a/a0/S6r3b8xsHNHfYAnwxUoNzR9aMRtc2sw+K6/nO2OOaXY3JGfOWXpNkHYaucyuHlSTUERyKcQc9Ozhx6Yac6Y+25HLj3ENXWZ399CjGnk6aQH7rryOF6dObHY3JGcGz74vSDtrg7TSPA0doPddeV0jTyctItSLUaTculzGxemlTfXeHLgU2IlozvkEomV2VVX1/tuOB2ToqhTRhxbdrtU90keoBQWtvooj7Troq4j2Ob3UzDYANgbOAF509/PMbAawhbt/q0I7moMWkVRCzEFfPeK4VGPOcV1X53IkT7PMblNgT+BzAO7+FvCWmR1KlPYIcBVwLzDgAH32NrpbL72duewantj+wGZ3Q3LmI0/dFqSddpji2I5ox7orzGwX4CHgVGqo6i2SJNSLUaRcqy+zSzNADwI+Dpzi7g+Y2fnAjLQnKN1n9Zld9tMUh/Sy3aN3K5NQ+giVvLS2xSPo9VIc0wl0uvsD8c+ziQbs7riaN2mreouINNK6lF95laYm4XNmttzMdnT3RUR1CJ+Iv6qq6r3do3dn7K4UkVK9pV7yPPimkXYd9ClAR7yCYzFwPFH0XVVVb32UlXJmxmfGHN7sbkjO/HrpjUHayb4OpLnS7mb3CDA+4VdVVfVWpCRJQr0YRcq1SwQdhCJoKaeq3pIk2E3CIK00j6p6S9PpupB6aYd10P2leu9PlVW9/zTisNp7KoU0oWsuiz+2b7O7ITkTakFBu0xxnA/c7u5TS1K996fKqt4TuubW0EUpOq3ukXop/AA9QKp31Sf7zfBjq36OFNunn+3QHLT0EWraq9WvrCyp3lBlVW+RJJqDlnpp9TnoirvZmdl44H5gj5JU71eIyotXVdWb1n9Dk8C0ikOSxNdF5uH1+2PS7WZ3+tJ87mZXc6q3u3e7+1p3XwdcAkyoVydFRGqxDk/1lVc1p3rXUtVbH2Ulia4LqZfC3ySMJaV6/1e1Vb31UVbKaYpDkugmYSRLqvenqz2ZIiVJoutC6qVdIuggFClJOUXQkiTUm/aaFq+yl2Yd9I5ExWF7bAecCfyaKovGKlKSJLoupF5ae3hOd5NwETAOwMzWB1YAc4iqqswvKRo7gwo1CRUpSTkz48kdVJNQevvwfweqSRikleapdopjH+Bpd19aS9FYRUqSJNSLUaRcyCV0ZnYA0bYX6wOXuvt5Zb9/L9HMwq7AKuBod1+S5ZzVDtDTgJnx96mKxpZSBC3lzIxNNt622d2QnFn92jNB2gk14sSzB78A9iPKDXnQzOa5+xMlh50IvOTu25vZNOAHwNFZzpt6gI6X2B0CnF7NCZTqLZWEejGKlAs4xTEBeMrdFwOY2SzgUKLSfz0OBc6Kv58NXGBm5hki02oi6AOBh929O/65uydZZaCisaVVvf+w9VSF0NLLJ7pv4IZh05vdDcmZI7pnVj4ohbXhpjhGAstLfu4E/rm/Y9x9jZm9DAwh2hKjJtUM0NN5d3oDYB5VFo39RPcNVXVO2kOoF6NIubQRdMIn/Yvj4PKdQxKeVj76pzmmKmk37N+YaO6lNFvwPKosGvvy5/erpY9SYJtdehdv/P0Pze6G5MyGH9wjSDuecnws/aTfj05gdMnPo4Cufo7pNLNBwGbAi6k7myBtJuFrRKF66WOrqLJorIhIIwWcg34Q2MHMtiVaajwNOKbsmJ5ZhT8CU4F7ssw/Q4MzCTe79K5Gnk5aRKhoSaRcqGV28ZzyycAdRMvsLnf3x83sbGChu88DLgN+Y2ZPEUXO07KeV6ne0lRK9ZYkedwsKa65emvZY2eWfP8GKaZ6q5El1Xtzqiwaq0QVSaLrQuplTYsne2dJ9T6eKovGKlKScoqgJUm4CLq1r60sqd5Vn0yRkiTRdSH10m57cZSmekOKorGllvyTFn1Ib2P/PF8RtPShCDqSJdX7QqJisT1FY/8T6FM0VqneUokiaKmXdoqge6V6l6R8Y2aXADcnPal0Afj/GXtUa7+dSXAXLr2e48cc0exuSM5csTRM1vHaFv90VnOqdy1FYy9cen11vZO2EOrFKFIuzxW708iS6v3DaovGfm+bY2voohTZGcs6NActfWgOOpIl1bvqorFnLOuo9inSBjQHLfXSTnPQmSlSknJmxupvTml2NyRnNvnhnCDttMUUh4hIK2qLKQ4z+yrweaL55seIsgiHA7OAwcDDwKfd/a0K7WTqrBRTqGhJpFzhV3GY2UjgX4GPuPvrZnYdUcLKQUSp3rPM7FdE9bguHKgtTXFIOaV6S5JQwVy7THEMAjYys7eBjYFngb15dz/Uq4hqcQ04QCuCliS6LqReCn+T0N1XmNmPiaqmvA7cCTwE/MPd18SHdRLV4+qjNJPw3JEfYvqQUSH6LQWx3aN3M32bw5rdDcmZmcvmBmmn8HPQZrYFUbXabYF/ANcTZRWWS/xLlGYSmpl/e8Xfau6sFFOoF6NIuXaY4tgXeMbdnwcwsxuBfwE2N7NBcRSdVJ+rjz9sfXiWvkoB7fHcjZqDlj6CJaq0+LWVZoBeBuwWZxO+TrTl6ELgd0R1t2aRsqr3Hs/dWHtPpbA0By31srboEbS7P2Bms4mW0q0B/kw0ZXELMMvM/iN+7LJKbV01Qqne0ttnuzrYdsi4ZndDcuaZVY8EaafVpziskR8BzKy1/1oi0jDunvmj1T6jPplqzJnfeWcuP8Yp1VuaysyYss2nmt0NyZk5y24K0k6rR9BK9RaRwir8MjvoN9X7V8BE4OX4sM+5+4ATR7oZJElCRUsi5do51RvgNHefnfZkj217cG29lMLa+ZlbeOXLBzS7G5Izm154e5B22mWKozzVu+Ka5yQ7P3NLLU+Tggv1YhQpV/gBOinV293vNLNjgHPN7ExgPjDD3d8sf35pqve3N/8Yh79/bMj+S4vbtXMeLxy0Z7O7ITmz5a0LgrTT6gsTKi6zi1O9bwCO5t1U79lEg/JzwAZE66KfdvezK7TV2n8tEWmYEMvsJoyYmGrM+VPXfbm8QVZzqre7Xx3//k0zuwL4RqWGWv3dTMIzM1YdpghaehsyN1AEXfQpDvpJ9e6p6m3R0ozDSFHVW6s4JEmoF6NIubXe2huOZkn1vs3MtgIMeAT4UqW2lOot5T7bpare0pc2S4oo1VtEcinEHPQuW/9LqjHnL8/9v1x+vFeqtzSVmWnDfumj1TbsN7MjiapKfRiY4O4L+zluCbAaWAuscffxA7Xb0AFac9CSRBv2S72sa1xQ+FfgcOCiFMfu5e4vpGk0bar3qcAXiOabL3H3n5nZYOBaYCywBDjK3V8aqJ0LR2kOWnr7cmcHp42ZVvlAaSs/WjorSDuNiqDd/UkIH4SuV+kAM9uJaHCeAOwCTDazHYAZwHx334E4USVoz0REMlrr61J9mdlJZraw5OukOnXJgTvN7KE050iTqHIksL+7fz7++TvAm8CJwKR4qd1w4F5337FCW5qEFpFUQtwk/OBW41ONOX9/fmHFc5nZ3cDWCb/6trv/Nj7mXuAbA8xBj3D3LjMbCtwFnOLu/a4zTTPF8VeilO4hROugDyIqeTXM3Z8FiAfpof106J1U74suuoiTTqrXG5O0IjPjZ6M19SW9/dvyjiDthJzicPd9A7TRFf+70szmEM1M9DtAp1pmZ2YnAl8BXgWeIBqoj3f3zUuOecndt6jQjiJoEUklRAT9gS0/nmrMefqFh4NMHg8UQZvZ+4D13H11/P1dwNnu3u9uYaluErr7ZcQ1B83se0An0F2STTgcWJminTSnkzZiZroupI9giSqNW2Y3Bfg5sBVwi5k94u77m9kI4FJ3PwgYBsyJ/2+DgGsGGpwhfQQ9NA7JtwHuBHYHzgBWuft5ZjYDGOzu36zQjl6JIpJKiAh6zJCPpRpzlq56NJdrgNOug74hnoN+G/iKu79kZucB18XTH8uAIys1okhJyimCliRK9Y6kneL43wmPrSLaOCk1JapIEl0XUi+F37A/pFZ/N5PwFEFLEkXQEaV6S9PpupB6aWCqd11kSfU+K37s+fiwM9z91oHaeWjUIRm6KkW0a+e8lo9yJLxWW8VRL2mqepemer8F3G5mPdVff+ruP057sl0759XUSSk2RdBSL4XfsJ9o+7z73f01ADO7D5hSy8kUKUk5zUFLEs1BRypulkSU6r2nmQ2Jy14dBIyOf3eymT1qZpfHxWX7KN2E5OKLLw7UbRGRyta5p/rKqyyp3ucBLxDtznQOMNzdT6jQTn7/EiKSKyESVbZ4//apxpyXXn0ql/NsVZe86kn1dvdfljw2FrjZ3Xeq8HQN0NKLpjgkSXxdZB40N3v/B1JdXC+/+nQuB+i0qzhKU70PB3bv2YcjPmQKquotNdJ1IfXS6m/+WVK9f2Nm44ii4iXAFys2Mmx6zR2VYjqieyZ/2/GAZndDcuZDiwbcQyi1Vl/FoareIpJLIaY4NtpoTKox5/XXl+byY1xDMwkXfWj/RvBS7gAAAAGeSURBVJ5OWsCOf7uDpbvu3exuSM6MeeieIO20+hSHImgRyaUQEfR7Nxydasx5843liqBfOXVyI08nLWDT82/m9QW/bnY3JGc22vMzQdpRBF3NyRRBi0hKISLoQRuMTDXmrHlrRS4j6IYO0PIuMzvJ3ZVaKb3oupBSaVK9pT5U3lyS6LqQd2iAFhHJKQ3QIiI5pQG6eTTPKEl0Xcg7dJNQRCSnFEGLiOSUBmgRkZzSAN1gZnaAmS0ys6fMbEaz+yP5EFclWmlmFbftlfahAbqBzGx94BfAgcBHgOlm9pHm9kpy4kpA+65KLxqgG2sC8JS7L3b3t4BZwKFN7pPkgLsvAF5sdj8kXzRAN9ZIYHnJz53xYyIifWiAbqykDVm0zlFEEmmAbqxOYHTJz6OArib1RURyTgN0Yz0I7GBm25rZBsA0YF6T+yQiOaUBuoHcfQ1wMnAH8CRwnbs/3txeSR6Y2Uzgj8COZtZpZic2u0/SfEr1FhHJKUXQIiI5pQFaRCSnNECLiOSUBmgRkZzSAC0iklMaoEVEckoDtIhITv0PQ/gdcPJVR7kAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Sample self-masked mnar data -----------------------------------------\n", "X_miss_selfmasked = produce_NA(X_complete_cont, p_miss=0.4, mecha=\"MNAR\", opt=\"selfmasked\")\n", "\n", "X_mnar_selfmasked = X_miss_selfmasked['X_incomp']\n", "R_mnar_selfmasked = X_miss_selfmasked['mask']\n", "\n", "print(\"Percentage of generated missing values: \", (R_mnar_selfmasked.sum()).numpy()/np.prod(R_mnar_selfmasked.size())*100, \" %\")\n", "\n", "ax = sns.heatmap(X_mnar_selfmasked, mask=R_mnar_selfmasked.numpy()==1, linewidths=0.005, linecolor='black')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### MNAR with quantile censorship" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With `opt = \"quantile\"`, missing not at random data are generated by using quantile censorship. First, a subset of variables which will have missing variables is randomly selected. Then, missing values are generated on the q-quantiles at random. Since missingness depends on quantile information, it depends on masked values, hence this is a MNAR mechanism.\n", " \n", "Two additional arguments are required: \n", "\n", "* `p_obs`: the proportion of fully observed variables. Note that at least one variable is chosen to be observed. \n", "\n", "* `q`: quantile level at which the cuts should occur." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Percentage of generated missing values: 16.5 %\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAD7CAYAAABHYA6MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de5xVdbnH8c+jZN7zkgzD1QtmmqYdeJFmJYl5JbwjGKVlmR2li5UhlbfSsJt21JO3TMsBRBQjCEXxlueoJzQ1L5moMAwDM96TNBV4zh9rYfs2s9fe+zd7r7339/168ZKZWfu3fvKa/cwzv/V7fo+5OyIikj4b1HoCIiJSmAK0iEhKKUCLiKSUArSISEopQIuIpJQCtIhISlUUoM3sYDN72syWmNmUUJMSERGwcvdBm9mGwN+BTwMdwJ+Bie7+ZLjpiYg0r34VvHYUsMTdnwMws5nA4UCPAdrMVBUjIom4u1U6xjsvPpco5rzn/TtWfK++UEmAHgQsz/i4A/ho7kVmdjJwMsBPd96Zzw0cWMEtpdEMuOce3rjxh7WehqTMpsf+oNZTSIVKAnShnzh5P63c/UrgSogy6O8880wFt5RGpDej9Jl1a2s9g4pUEqA7gCEZHw8GOnt7wetTjqrgdtKItph2M6tGf7LW05CUGXD3vWEGWrsmzDg1UslDwn5EDwnHACuIHhIe7+5P9PIarUGLSCIh1qDf7nwiUczZaOCHGmsN2t3XmNlpwG3AhsA1vQVngPtaji73dtKgPt51EzpRUXKZBYqX69aFGadGys6gy7qZMmgRSShIBr380WQZ9JA9GyuDLscDrUdW83ZSB/ZeOUcZtOQJl0E370NCEZF08/pe4qgoQJvZUuB1YC2wxt1H9nb93ivnVHI7aVDBsiWRHF7nuzhCZNCfcvcXk1x4T8sxAW4njWS/rtlc3/rZWk9DUmbSyrYwAzXzQ8I4gx6ZNEDrIaGIJBXiIeFbf78vUcx57wc+nspf4yrNoB1YGAfeK+KqwSyZpd4X7vgBJrWo1Fv+bdD9d3NLy8RaT0NS5oiuGWEGqvOHhJVm0APdvdPM+gO3A5PdvccSIGXQIpJUkAz6qbuSZdC7fqrxMmh374z/221mc4hOuOsxQD80eFwlt5MGNKJjrrbZSZ5gD46b9SGhmW0GbODur8d/PxA4r7fXjOiYW+7tpIFpF4f0mTp/SFhJBt0CzInfXP2A6e5+a28vUAYtuUZ0zOXJ4YfUehqSMrstWRBkHPcmXoMu+WZagxaRhEKsQf/rkXmJYs7Ge41N5a9xVa0k/OsOh1XzdlIH9nh+vtagJY8OS4pUNUDv8fz8at5O6oTWoKXPNHqpt5ldA4wFut199/hz2wA3ANsDS4Hx7v5KsbGuHKSKMcl28oo2ntpZa9CSbddnwqxBs/adMOPUyAYJrrkWODjnc1OARe6+M7Ao/lhEJF3WrUv2J6USPSQ0s+2BeRkZ9NPAaHdfaWatwN3uvkuCcbTYKCKJBHlIeP+MZA8J95mYynW2ctegW9x9JUAcpPv3dGFmqfdnthnFyM2Hl3lLaURnt0/XQ0LJo4eEkT5/SJjb1fsPL/9fX99S6oweEkqfadIA3WVmrRlLHN1JXqRDcSTXEV0zePXz+9d6GpIyW/32ziDjeJ0/JCw3QM8FTgCmxf/9fZIXBTuhShpKqDejSJ4m2GY3AxgNvN/MOoCziQLzLDM7CWgHjk1ys0X9E10mTWRM9428/cKztZ6GpMxG2+0UZqA6X+JQqbeIpFKIXRxv3nF5opizyQGnpPJBSFUrCW/SGrTkOLprBgv7j6/1NCRlDuyeFWagOs+gqxqgj9YatBQQ7M0okqsJ1qALlXqfA3wZeCG+bKq7/7HYWP/4ykHlz1Qa0pZX3MY/p51Y62lIymw25dowA60Je2C/mW0ILAZWuPvYnK+dCPwUWBF/6lJ3v7qS+yXJoK8FLgV+m/P5i9z9Z6XcbMsrbivlcmkSwd6MIrnCZ9BfB54Ctuzh6ze4+2mhblY0QLv7vXGpd8X+Z8BRIYaRBrLvqpu1u0fyjOm+McxAAdegzWwwcBhwPnB6sIF7keSwpJ6cZmaPmdk1ZrZ1TxeZ2clmttjMFv/+jecruJ2ISIl8XaI/mXEq/nNygdEuBs4Aeov6R8dxcbaZDal0+uUeltQCvAg48EOg1d2/mGAcbbMTkUSCbLObMy3ZNrsjp/R6LzMbCxzq7v9pZqOBbxdYg94WWO3ub5nZKUTHMFdUJlvWLg5378qY1FXAvCSvGztEHVUk27zl8/nWsAm1noakzM+XzQwzULg16H2BcWZ2KLAxsKWZXe/uk969lftLGddfBVxY6U3LCtDrz+GIPzwSeDzJ6+YtV0cVyRfszSiSK9AuDnc/EzgTICODnpR5TU5cHEf0MLEi5ZZ6jzazvYiWOJYCX0lyszkqVJEcR3bN4BdD1GlHsp2+vC3MQH1cKW1m5wGL3X0u8DUzGwesAV4GTqx4fJV6i0gaBVmDnnF2sjXoieeq1Pvlo/ar5u2kDmxz8z06sF/y6MD+SFUD9DY331PN20md0IH90meaoNR7CFEV4QCi/X9Xuvsvy+nsfbW6ekuOL61o42KtQUuOb4Rag167Nsw4NVJ0DTrumNLq7g+b2RbAQ8ARRAvgL7v7NDObAmzt7t8tMpZ+lxWRRIKsQf/mjGRr0F/4SSp/jUtS6r0SWN8g9nUzewoYBBxOtLsD4DrgbqDXAD1p6JEVTFUa0fXtc7QGLXm0Bh0pqdQ7rij8CPAgOZ29gYKdvTNLKJ9ZrVJvEamihKXeaZX4IaGZbQ7cBHzD3f+R9CdcblfvB19+pJx5SgPTQ0LpK76uvn87SxSgzew9RMG5zd1vjj9dcmfvhwaPK3+m0pBGdMxlw/cMrPU0JGXWvtMZZqA6X+JIsovDgF8DT7n7LzK+VHJn7xEdc8ucpjSyYG9GkVx1vosjSQa9L/A54K9mtn59YipldPa+ZLC2U0m2yR1tXDBU3xeSbWp7oG12dZ5Bq9RbRFIpxDa7N355SqKYs+nXL0/lgxB19ZaaOrprBiMGfLzW05CUeWjVfWEGqvMtnOrqLTUX7M0okqvOlzgqKfU+hxI7e08dpgxasl2wbIYKVSRPuEKV+v7eSpJBrwG+lVnqbWa3x18rqbP3BcuUQUs+7YOWPtPouzh6KfUu2cwBx5fzMmlgE1ZNVwYteUL90PZGX+LIlFPqvS9RZ+/PA4uJsuy80+zi7riFOuSKAMqgpQ81wRIHULDU+1dEHb3Xd/b+OZDX2Tuz1PuM7SfW97+WBPfTZTPV8kryhGt51QQZdKFS73I7e4uIVE2jZ9A9lXqX09n7p+reLAUEy5ZEcq1p8IeE9FzqPbHUzt73thxd5jSlUX2y6yZu1MNjyXHsqulhBqrzJQ6VeotIKoUo9f7n945NFHM2O//GVD6prmol4T0tx1TzdlIH9uuazYL+x9V6GpIyh3TfEGScht9mZ2YbA/cC742vn+3uZ5vZDsBMYBvgYeBz7v52b2Pt1zW78hlLwwn1ZhTJ0+gPCYG3gP3dfXW8m+M+M1sAnE5USTjTzC4HTgJ+1dtA/6XjRiXH1zramKIjACTHtFBVx3UeoEtagzazTYH7gK8C84EB7r7GzPYBznH3g4q8vr7/tUSkakKsQa8+fVyimLP5L+bW7xq0mW0IPAQMBy4DngVedfc18SUdJCj/VgYtub7W0cZ3hk2o9TQkZUJtyW2KnoTuvhbYy8y2AuYAuxa6rNBrVeotxWh/vPSZZgjQ67n7q2Z2N7A3sJWZ9Yuz6MFAwcZymaXeE4cdUd//WhLczPbfaw1a8oRbgw6zi6OnzRI517yX6GjmEcBLwHHuvrSS+ybZxbEd8E4cnDcBDgAuBO4CjiHayZGoaezM9qKXSBMK9mYUyRUugy64WcLdH8i45iTgFXcfbmYTiOJkRXtIk2TQrcB18Tr0BsAsd59nZk8CM83sR8BfiMrBezVLFWOSY/yq6TzQemStpyEps/fKOWEGChSgPdpNsTr+8D3xn9zBDwfOif8+G7jUzMwrqAZMch70Y0RHjOZ+/jlgVLk3FhHpa742XKFK7mYJd38w55JBwHKAeHfba8C2wIvl3rOqlYTjQ9XXS0MJli2J5EqYQRfYzHBl/PzsXbmbJcxsd3fPPCSu0Fa9ilL4qgboy7TNTnKc2tFGy/s+WOtpSMp0vfa3IOMk3WaXuZkhwbXrN0scTPYpnh3AEKDDzPoB7wNeLmW+uSop9b4W2A94Lb70RHd/pPAokVM7dKyk5Av1ZhTJE2gNupfNEpnmEm2YuJ9oA8Wdlaw/Q2Wl3gDfcffEB2xcOUgZtGQ7eUUbnx2qh4SSra091EPCMMPQ82aJ84DF7j6XaKPE78xsCVHmXHEFVpKHhEmeXiZy8gpl0JIv2JtRJIevCROhe9kscVbG3/8FHBvkhrGySr3d/UEz+ypwvpmdBSwCprj7W72N89DgcZXOVxrMiI656uoteYI1Eq7v00bLK/U2s92BM4FVwEZEi+vfBc7Lfa1KvaUYdfWWvtIUZ3Gsl/n00t1/Fn/6LTP7DfDtHl7z7tPRtoGT6vtfS4KbtLKN1eeqgEmybX52oC25jZ5B9/T0cn3T2Lip7BEkaBo7aaXWoCVfsDejSI5myKB7enp5Zxy8DXgEOKXYQNe3aheHZJu0so0lHzqw1tOQlBn+xMIwA9V5Bq2msSKSSiEO7H/psP0SxZxt59+TygchVa0kPG+o1hol21nt0/nSMDUTlmxXLwvTv9TrPIOuaoAWEamqOg/QiZc44jXoxcAKdx9bTldvLXGISFIhljhe+HSyJY7tbq//JY6vA08BW8YfX0iJXb1valHnDMl2dNcMLtEhWpJjcqBze+p9iSNRBm1mg4HrgPOB04HPAC+grt4i0kdCZNBdo0cnijktd99d1xn0xcAZwBbxx9uSsKt3ZiXhaVuM5OBNhpc/W2k4Y7tnqtRb8oSqLq33DDpJocpYoNvdHzKz0es/XeDSgu+yzEpCM/NLX19c5lSlUanUW/qKr6vv760kGfS+wDgzOxTYmGgN+mISdvXO9CNts5Mc32+fTue+o2s9DUmZgf9zd5Bx6j2DLqlQJc6gvx3v4rgRuCnjIeFj7v7fRV6v32VFJJEQa9Ar9tk/UcwZdP+dqUy1K9kH/V1K7Oq9eqoKEiTb5hfM1hq05NEadESl3iKSSiEy6PaRYxLFnKGLFzVcBl0yHdgvuXRgvxQSLoNOZdxNTKXeItKwmiZAFyj1vpYSu3qP6Jhb7jylgWmbnfSVev/lrJJSbyixq/ei/kH7KUoDGNN9o5Y4JI+WOCJJm8YOBg7j36XeZRnTfWO5L5UGpgxa+kqA54w1VW6p93pFu3pnlnqfte0eHLvl0AqmK41mj+fn84NhKmCSbD9cFqYN2tq19R2gi26zi0u9D3X3/8wpVGklu6v3s+6e19U7Zyz9LisiiYTYZvf0Bw9JFHN2+duCVEbyskq9zex6d58Uf73Xrt6Znv5gr4fdSRPa5W+3cUf/8bWehqTMAd2zgoxT72vQlZR6Z3b1vgj4l7tPKfJ6ZdAikkiIDPqpnQ9NFHN2feaPqYzkleyDbiu1q/eTww+p4HbSiHZbskAFTJIn1JbcpsqgK76ZMmgRSShEBv3XHT6TKObs8fwfUhnJq1pJ2LX/J6t5O6kDLXfey28GquWVZPtCZ6CWV3WeEiqDFpFUCpFBPzJsXKKYs9eyuUXvZWbXAOsbmOxe4Oujgd8Dz8efurnYzrZikhaqLAVeB9YCa9x9pJltA9wAbA8sBca7+yu9jXN/65GVzFUa0D4r57BQuzgkx4GhdnGELVS5FrgU+G0v1/zJ3ceGuuEGJVz7KXffy91Hxh9PARa5+87EhSqhJiUiEoJ7sj/JxvJ7gZf7dMI5KlmDPhwYHf/9OuBuokP8e7TPyjkV3E4aVahsSSTXuoQZdGbFc+zKuJ9qqfYxs0eJWgB+292fKGOMdyUN0A4sjNeQr4gn3uLuKwHi/dD9C70wt6v3IZvsVMl8pcEc1n0Dbz6oAC3ZNvlomGWvteuSLRJkNreuwMPAMHdfHRf23QLsXMmAiR4SmtlAd++Mg/DtwGRgrrtvlXHNK+6+dZFx9JBQRBIJ8ZDwgYFHJYo5e3fenOheZrY9MK/QQ8IC1y4FRrr7i0nGLiRRBu3unfF/u81sDjAK6MqoJmwFuouNM6//hHLnKQ1qbPdMHTcqeUKdcJh0iSMEMxsAdLm7m9koomd8L1UyZtEAbWabARu4++vx3w8EzgPmAicA0+L//r7YWGO7Z1YyV2lQOm5U+krIXRxmNoPoudv7zawDOBt4T3Qfvxw4Bviqma0B3gQmeIXZR5IMugWYE7+J+gHT3f1WM/szMMvMTgLagaKn8f9tl4Mrmas0oA8+fasyaMkTLIMOMkrE3ScW+fqlRNvwgikaoN39OWDPAp9/CRhTys0++PStpVwuTUIZtPQVp76/t6pa6n3hUJX0SrbvtrdxyWB9X0i2yR1hSr3X1HlHFZV6i0gqhdjFsajluEQxZ0zXDamM5JWUep8DfBl4Ib5sqrv/sbdxXj9Dpd6SbYufzOGNm86v9TQkZTY9+ntBxgm5Bl0LSfdBLyVnP18coFe7+88S30wZtIgkFCKDXtgyIVHMObBrZv1m0KG8cvynqnk7qQNbT7+LN27+ca2nISmz6VFnBhmn3jPopIclrS/1figu3V7vNDN7zMyuMbOCVYRmdrKZLTazxdcu6ax4wiIiSa3FEv1Jq0pKvZ8GXiQK3j8EWt39i0XG0RKHiCQSYonjDwMmJoo5n1k1I5VRuuxS7/joPQDM7CpgXrFxnt390+XOUxrUTo/frkIVyROuUCWVcTexsku915/DEV92JPB4sbF2evz2iiYrjUmFKtJX6v1HfyWl3r8zs72I/g2WAl8pNtD/DTyigqlKIxrVeQsXD1GhimT7xvIwhSr1/pBQhSoikkoh1qBnt342Ucw5ZmVbKn+Nq+o2uwfUk1By7L1yjtagJU+oZa+1QUapnaoG6L3V8koK0Bq09JV1df6tlbTUeyvgamB3ojXnLxJtsyupq/fSj5R0+J00ge3/ski7eyRPqA0F9b6LI+k+6OuI2olfbWYbAZsCU4GX3X2amU0Btnb3XpvGag1aRJIKsQZ9/cBJiWLOpM7rUxnJk2yz2xL4JHAigLu/DbxtZiV39daxkpJrckcbz334gFpPQ1Jmx8fuCDJOMyxx7Eh0Yt1vzGxP4CHg65TR1VukkFBvRpFc9b7NLkmA7gf8BzDZ3R80s18CU5LeILOd+Yp99tcSh2QZ/MBd3NNyTK2nISmzX9fsIOOsrfMMOslhSR1Ah7s/GH88myhgd8XdvEna1VtEpJrWJfyTVkl6Eq4ys+Vmtou7P03Uh/DJ+E9JXb0HP3BXhdOVRhQqWxLJlebgm0TSfdCTgbZ4B8dzwBeIsu+SunpfPkgPCSXbKSvaVKgieULtja/zloSJT7N7BBhZ4EslbWw+ZUWY+nppLCpUkb7SLBl0EB17q6OKZBv8wF18c9hxtZ6GpMxFy24IMo5KvUugNWgpJNSbUSRXM+yD7qnU+yBK7Or91x0OK3+m0pD2eH6+frOSPKGSuWZZ4vglcKu7H5NR6n0QcFEpXb33eH5+GVOURqffrKSvNHyA7qXUu+SbzW2ZUPJrpLGN65qpXRySJ9gujiCj1E4lpd4QdfX+PLAY+Fah0+xU6i3FaBeH9JWQa9BmdjDRasKGwNXuPi3n6+8FfguMAF4CjnP3pZXcs5JS70uJunmv7+r9c6K16SyZpd6/HZTsZClpHid0tvEN7eKQHBenbBeHmW0IXAZ8mqi6+s9mNtfdn8y47CTgFXcfbmYTgAuBir65yy71dvcud1/r7uuAq4BRlUxERCS0dXiiPwmMApa4+3PxMu9M4PCcaw4nOtkTojg5xir89bDsUu9yunqf0KlCFckXKlsSyZX0IWGBpdgr49/+1xsELM/4uAP4aM4w717j7mvM7DVgW+DFkiadoZJS7/8qtav3wv7jy5ymNKoDu2cxeZi+LyTbJctmBRkn6Zpq5lJsDwplwrnDJ7mmJJWUen+u1Jsd2B3mH10aS6g3o0iugNvsOoAhGR8PBjp7uKbDzPoB7wNeruSmVa0kvEMZtOQ4oHuWttlJnlA7e9aE67L3Z2BnM9sBWAFMAI7PuWYu0cme9wPHAHd6hd/cSfZB70LUHHa9HYGziLaTlNQ09gBl0FKAttlJXwkVnuM15dOA24i22V3j7k+Y2XnAYnefC/wa+J2ZLSHKnCsu/EjUNPbdi6OtJiuIFsdPpcSmsZcO0TY7yTa5o43n91RPQsm2w6N3BGkae+b2xyeKOT9eOj2VWUKpSxxjgGfdfVk5TWMnd2gXh+Tb4VH1JJS+kXALXWqVGqAnADPivydqGptJB/ZLLh3YL4Wo1DuSOEDHW+zGAWeWcgOVeksxWoOWvtLwhyVlOAR42N274o+71her9NY0NnN/4SPDxtX7DzQJ7CPtf9DuHskTakPB2jrPoUsJ0BP59/IG/HtLSeKmsR9p/0NJk5PmoN090leaIoM2s02JDgnJrBacRolNY/8xWQf2S7YtL5nPv5bcX+tpSMpsPHyfION4M2TQ7v4GUU155udeosSmsSIi1dQUGXQoW16ijiqSL1S2JJKr2bbZVeS+lqOreTupAx/vuolxQ7T0JdnmLg+TzNV3eK6s1HsrSmwa+/Gum8qcpjSyUG9GkVxr6jxEJzkP+mlgL8gq9Z5DdORoSU1jpw6bWOY0pVFdsGwGHxuort6S7X87wzQSboqHhBkyS71LvtkFy2YUv0iaTqg3o0iuZntImFnqDQmaxmZa+Yn9SrydNLrWP92jIwAkzykrwpzbU+8ZdOLT7OJS707gQ+7eZWYtRK1c1jeNbXX3vKaxOaXeI4LMWkQaXojT7E7Y/uhEAe66pTel8ryBsku9M0q+MbOrgHmFXpRZ6v39hEf/SfM4f9kMvjWs4mNzpcH8fNnMIOOsrfODuMou9S6naez5WoOWAkK9GUVyNcU+6B5KvX9SatPYXw3WWqNk+2pHG3NblEFLtnFdYX5oN80adJCbWbgGYSLS2EKsQR837IhEMeeGZbfU/Rp0xWYNyO2xKM1u/KrpOrBf8oQ6I7wpljhEROpRvS9xJF2D/ibwJaL15r8SVRG2AjOBbYCHgc+5+9u9jTN+1fSKJiuNSR1VpK80/C4OMxsEfA3Yzd3fNLNZRAUrhxKVes80s8uBk4Bf9TbW2UO1xCHZzm2fziFDDqn1NCRlFixfEGScZlni6AdsYmbvAJsCK4H9gfUR9zrgHIoE6HPblUFLvlBvRpFcDV/q7e4rzOxnRF1T3gQWAg8Br7r7mviyDmBQoddnVhJO22EXJrUMDDFvaRCDH7iLU4epJ6Fku2xZmDZoDb8GbWZbA4cDOwCvAjcSVRXmKvgvkVlJaGY+5fmny56sNKZQb0aRXM2wxHEA8Ly7vwBgZjcDHwO2MrN+cRY9mOicjl49OmxsJXOVBrTnsnm8Omn/Wk9DUmar6+8MMk69b+FMEqDbgb3jasI3iY4cXQzcBRxDtJMjUVfvPZcVPK5DmlyoN6NIrrWNnkG7+4NmNptoK90a4C9ESxbzgZlm9qP4c78uNtYtLTqwX7Id0TWDka2fqPU0JGUWr/xTkHHqfYlDpd4ikkohSr3HDD4wUcxZ1LEwlZvxq1pJOHHoEdW8ndSBGe231P06oYRXb6XeZnYs0VbjXYFR7r64h+uWAq8Da4E17j6yt3FV6i0iDauK2+weB44Crkhw7afc/cUkg1ZS6n05sB/wWnzZie7+SG/jzGi/JcntpMmo1Fv6SrVKvd39KQj/vVxJqTfAd9x9dtKbLfnQgeXNUhrW8CcWaolD8lR7iSOnNR/AlXENR2gOLIyfx11R7B7llnoX3fNcyPAnFpbzMmlwyqClryQN0JkFdT0xszuAAQW+9D13L7rNOLavu3eaWX/gdjP7m7vf29PFZZV6u/tCMzseON/MzgIWAVPc/a0C/1Pv/mQ6t/+HGP++oQn/P6QZ7PrMAmXQkifUD+2Q31vufkCAMTrj/3ab2RxgFFB+gC5U6m1mk4AzgVXARkQ/eb4LnFdgQlml3md3P1Hi/5I0OmXQ0lfStA/azDYDNnD31+O/H0iBmJmp7FJvd78+/vpbZvYb4NvFBprXX73nJNvY7pnKoCVPsAy6etvsjgQuAbYD5pvZI+5+kJkNBK5290OBFmBO/P/WD5ju7rf2Nm7Zpd7ru3pbdLcjSNDVe2y3ujdLPmXQ0lfWenUOHHX3OcCcAp/vJDo7H3d/DtizlHErKfVeYGbbAQY8ApxSbCyVekuuI7pmKIOWPGlcg64FlXqLSCqFKPXec8DHEsWcR1f9byp/jatqJaFaXkmuc9vV1Vvy1dsadF+paoBWyyspRGvQ0lfW1fkP/6Sl3l8Hvky03nyVu19sZtsANwDbA0uB8e7+Sm/jTG9VBi3Zjl85nQuGfrbW05CUmdreFmSces+gNyh2gZntThScRxE9gRxrZjsDU4BF7r4zcaFKX05URKRUa31doj9pVfQhYXyM3kHu/qX44x8AbwEnAaPjrXatwN3uvkuRser7x5mIVE2Ih4Qf2G5kopjz9xcWp3KdLckSx+NEJd3bEu2DPpSo5VWLu68EiIN0/0Ivziz1HrPNSD68xU5BJi6N4aJlN3DdQC1xSLYTOrXEAQm32ZnZScCpwGrgSaJA/QV33yrjmlfcfesi49T3v5aIVE2IDHqn9/9Hopjz7IsP120Gjbv/mrjnoJldAHQAXRnVhK1Ad7FxDhlySCVzlQa0YPkCfjFEGbRkO325MmhInkH3j09fGgosBPYBpgIvufs0M5sCbOPuZxQZp77/tUSkakJk0MO2/XCimLPspcfqN4MGborXoN8BTnX3V8xsGjArXv5oB44tNoi2U0muqe1tKlSRPCr1jqjUW0RSKUQGPXib3RPFnI6XH6/rDDqIywYrg5Zsp3a0cfDgg2s9DUmZWzt6PYUzMWXQpdxMGbSIJBQig27dardEMWflq0/WbwbdQ6n3OfHnXogvm+ruf3xkyfwAAAL4SURBVOxtnCeHaxeHZNttiVpeST4dlhRJ0vIqs9T7beBWM5sff/kid/9Z0pvttmRBWZOUxqbDkqSvpLmMO4kkGfSuwAPu/gaAmd0DHFnOzX6lNWjJ8dUO7eKQfNrFESl6WBJRqfcnzWzbuO3VocCQ+GunmdljZnZN3Fw2j5mdbGaLzWzxn1Y/E2jaIiLFrXNP9CetKin1nga8CDjwQ6DV3b9YZJz0/kuISKqEeEi49ebDE8WcV1YvSeU6W8m7ONaXerv7f2d8bntgnrvv3ttrzx/2WQVoyfL99uks6l+0xkmazJjuG4ME6PdtvlOimPPa6mfrN0D3UOq98frT7Mzsm8BH3X1CkXEUoEUkkRABesvNdkwUc/7xz+dSGaArKfX+nZntRbTEsRT4SrFB7ug/vuyJSmM6oHsWSz8yptbTkJTZ/i+LgoxT77s4VKgiIqkUIoPeZJNhiWLOm28uq+sMOohlI/av5u2kDgx76E5WfmK/Wk9DUqb1T/cEGafet9kpgxaRVAqRQb934yGJYs5b/1quDPr1M8qqb5EGtsVP5tR9liPhqVAlogxaRFIpRAbdb6NBiWLOmrdXpDKDrmqAln8zs5Pd/cpaz0PSRd8XkilJqbf0jZNrPQFJJX1fyLsUoEVEUkoBWkQkpRSga0frjFKIvi/kXXpIKCKSUsqgRURSSgFaRCSlFKCrzMwONrOnzWyJmU2p9XwkHeKuRN1m9nit5yLpoQBdRWa2IXAZcAiwGzDRzHar7awkJa4FDq71JCRdFKCraxSwxN2fc/e3gZnA4TWek6SAu98LvFzreUi6KEBX1yBgecbHHfHnRETyKEBXV6EDWbTPUUQKUoCurg5gSMbHg4HOGs1FRFJOAbq6/gzsbGY7mNlGwARgbo3nJCIppQBdRe6+BjgNuA14Cpjl7k/UdlaSBmY2A7gf2MXMOszspFrPSWpPpd4iIimlDFpEJKUUoEVEUkoBWkQkpRSgRURSSgFaRCSlFKBFRFJKAVpEJKX+H/8on+jEbAQFAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Sample mnar missing data with quantiles -----------------------------------------\n", "X_miss_quant = produce_NA(X_complete_cont, p_miss=0.4, mecha=\"MNAR\", opt=\"quantile\", p_obs=0.5, q=0.3)\n", "\n", "X_mnar_quant = X_miss_quant['X_incomp']\n", "R_mnar_quant = X_miss_quant['mask']\n", "\n", "print(\"Percentage of generated missing values: \", (R_mnar_quant.sum()).numpy()/np.prod(R_mnar_quant.size())*100, \" %\")\n", "\n", "ax = sns.heatmap(X_mnar_quant, mask=R_mnar_quant.numpy()==1, linewidths=0.005, linecolor='black')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "257.1875px" }, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 4 }