Multilevel models in Python, R and Mplus

bp_data.csv

mplus_bp_data.csv

Multilevel_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multilevel Modeling Software Comparison\n",
    "\n",
    "----\n",
    "David C. Folch (@dfolch)\n",
    "\n",
    "Reilley Luedde (@rluedde)\n",
    "\n",
    "\n",
    "## Background\n",
    "\n",
    "`R`, `Python`, and `MPlus` have multiple options for creating multilevel models. There is overlap among the models that can be built with these software, but not all models can be estimated in all software. In the following sections, you will see what needs to be done to get relatively congruent models.\n",
    "\n",
    "----\n",
    "## Software and Packages\n",
    "\n",
    "| Software | Package/Library |                                                    Class                                                    |\n",
    "|----------|:---------------:|:-----------------------------------------------------------------------------------------------------------:|\n",
    "| Python   |   statsmodels   | [MixedLM](https://www.statsmodels.org/dev/generated/statsmodels.regression.mixed_linear_model.MixedLM.html) |\n",
    "| R        |       lme4      |               [lmer](https://www.rdocumentation.org/packages/lme4/versions/1.1-23/topics/lmer)              |\n",
    "| Mplus*   |       N/A*      |                                                     N/A*                                                    |A                                                     |         \n",
    "\n",
    "*Multilevel modeling is available in the normal and demo version of Mplus. There's no need to install anything else.\n",
    "\n",
    "Both `Python` and `R` require the user to install external packages, statsmodels and lme4, respectively.\n",
    "\n",
    "## Estimators\n",
    "\n",
    "Each of the three software have various algorithms that are used to estimate the paramaters and standard errors. Below is a table of the estimation algorithms that we comparing. \n",
    "\n",
    "| Software | Available Estimators                                                                |\n",
    "|----------|-------------------------------------------------------------------------------------|\n",
    "| MixedLM  | Restricted Maximum Likelihood (REML)*, Maximum Likelihood (ML)                      |\n",
    "| lmer     | Restricted Maximum Likelihood (REML)*, Maximum Likelihood (ML)                      |\n",
    "| Mplus    | [See table in ESTIMATOR section](https://www.statmodel.com/HTML_UG/chapter16V8.htm) |\n",
    "\n",
    "                                   * = Default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What estimators produce similar models? \n",
    "\n",
    "By default, `lmer` (`R`) and `MixedLM` (`Python`) both use `REML`, and are accordingly very similar. However, there are **slight** differences in the calculation of the standard errors/deviations of these models. `REML` isn't possible in `Mplus` for TwoLevel random models. Maximum Likelihood (`ML`) is available in `R`, `Python`, and `Mplus` and the model results from using this estimator are very similar between all three software."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R setup for formatting output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compare each software, we will need to run `R` in this Jupyter Notebook. The `rpy2` extension allows this.  If you don't have it installed a `pip install rpy2` or `conda install rpy2` should do the trick.\n",
    "\n",
    "Once we load the `rpy2` extension, any cell that starts with the magic function `%%R` will be able to run `R` code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rpy2.ipython extension is already loaded. To reload it, use:\n",
      "  %reload_ext rpy2.ipython\n"
     ]
    }
   ],
   "source": [
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use:\n",
    "- R's multi-level modelling package `lme4`\n",
    "- `tidyverse`\n",
    "- `MplusAutomation` will help us extract parameters from Mplus for comparison with R and Python\n",
    "\n",
    "`install.packages(\"package\")` for any of these three packages that you don't have."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need some functions to extract paramaters from the models once they've all been built:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# Get a df containing the parameters for the fixed effects of the model.\n",
    "fixed_fx = function(R_model){\n",
    "\n",
    "    sum = summary(R_model)\n",
    "\n",
    "    fixed_fx = sum$coefficients %>% \n",
    "        data.frame() %>% # convert to df\n",
    "        select(-3) %>% # remove t-values\n",
    "        gather %>%\n",
    "        select(-1)\n",
    "\n",
    "    rownames(fixed_fx) = c(\"Int\", \"BMXBMI\", \"female\", \"Int_SE\", \"BMXBMI_SE\", \"female_SE\")\n",
    "    colnames(fixed_fx) = \"lme4\"\n",
    "\n",
    "    return (fixed_fx)\n",
    "}\n",
    "\n",
    "# Get a df containing the parameters for the random effects of the model.\n",
    "random_fx = function(R_model){\n",
    "    \n",
    "    sum = summary(R_model)\n",
    "    random_fx = sum$varcor %>%\n",
    "        data.frame() %>%\n",
    "        gather %>%\n",
    "        select(-1)\n",
    "    \n",
    "    random_fx = data.frame(random_fx[7:10,])\n",
    "    rownames(random_fx) = c(\"SEQN_Var\",\"Resid_Var\",\"SEQN_SD\",\"Resid_SD\")\n",
    "    colnames(random_fx) = \"lme4\"\n",
    "    return (random_fx)\n",
    "}\n",
    "\n",
    "lme4_params = function(model){\n",
    "    fixed_fx = fixed_fx(model)\n",
    "    random_fx = random_fx(model)\n",
    "    all_params = rbind(fixed_fx, random_fx)\n",
    "    write.csv(all_params, \"r_df.csv\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function reads in the output generated by `MPlus` and extracts certain paramaters to compare with `R` and `Python`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "mplus_params = function(mp_out_fp) {\n",
    "  mp_model = readModels(mp_out_fp)\n",
    "  coef_df = coef(mp_model)\n",
    "  coef_df = coef_df %>% \n",
    "    select(c(\"est\", \"se\")) %>% \n",
    "    gather() %>% \n",
    "    select(c(\"value\")) \n",
    "  \n",
    "  order = c(3,1,2,8,6,7,5,4,10,9)\n",
    "  coef_df = coef_df[order,,drop = FALSE]\n",
    "  colnames(coef_df) = c(\"MPlus\")\n",
    "  \n",
    "  rownames(coef_df) = c(\"Int\", \"BMXBMI\", \"female\", \"Int_SE\", \"BMXBMI_SE\", \"female_SE\",\n",
    "                       \"SEQN_Var\",\"Resid_Var\",\"SEQN_SE\",\"Resid_SE\")\n",
    "  write.csv(coef_df, \"mplus_df.csv\")\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python setup "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* If you have problems with the notebook, try `pip install html5lib==0.9999999`, because some versions of `html5lib` have a bug."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def sm_params(results_obj):\n",
    "    row_labels = [\"Int\", \"BMXBMI\", \"female\", \"Int_SE\", \"BMXBMI_SE\", \"female_SE\", \"SEQN_Var\", \"Resid_Var\", \"SEQN_SE\",\n",
    "                 \"Resid_SD\"]\n",
    "    column = [\"statsmodels\"]\n",
    "\n",
    "    results = pd.DataFrame(index = row_labels, columns = column)\n",
    "    \n",
    "    # HTML Parsing that might break depending on your version of html5lib\n",
    "    # Read results_obj into a list of DataFrames and select the first df.\n",
    "    summary_df = pd.read_html(results_obj.summary().as_html(), header = 0)[1]\n",
    "    summary_df.index = summary_df.iloc[:,0]\n",
    "    summary_df = summary_df.drop(summary_df.columns[0], axis=1)\n",
    "    \n",
    "    # Extract parameters \n",
    "    results.loc[\"Int\"] = summary_df.loc[\"Intercept\", \"Coef.\"]\n",
    "    results.loc[\"BMXBMI\"] = summary_df.loc[\"BMXBMI\", \"Coef.\"]\n",
    "    results.loc[\"female\"] = summary_df.loc[\"female\", \"Coef.\"]\n",
    "    results.loc[\"Int_SE\"] = summary_df.loc[\"Intercept\", \"Std.Err.\"]\n",
    "    results.loc[\"BMXBMI_SE\"] = summary_df.loc[\"BMXBMI\", \"Std.Err.\"]\n",
    "    results.loc[\"female_SE\"] = summary_df.loc[\"female\", \"Std.Err.\"]\n",
    "    results.loc[\"SEQN_Var\"] = summary_df.loc[\"SEQN Var\", \"Coef.\"]\n",
    "    results.loc[\"Resid_Var\"] = results_obj.scale\n",
    "    results.loc[\"SEQN_SE\"] = summary_df.loc[\"SEQN Var\", \"Std.Err.\"]\n",
    "    # results.loc[\"Resid_SD\"] = This isn't even in the normal printout ...\n",
    "    return results\n",
    "\n",
    "\n",
    "# Output a df that has selected lme4 results next to the statsmodels results\n",
    "def compare_params(r_csv_fp, results_obj, mp_csv_fp = None):\n",
    "    sm_results = sm_params(results_obj)\n",
    "    lme4_results = pd.read_csv(r_csv_fp, index_col = 0)\n",
    "    compare_results = pd.concat([lme4_results, sm_results], join=\"outer\", axis = 1)\n",
    "    # If we want to compare R, Python, and MPlus\n",
    "    # Since MPlus doesn't have a REML option, it's useful to compare only the Python and R REML\n",
    "    # models.\n",
    "    if mp_csv_fp:\n",
    "        mp_results = pd.read_csv(mp_csv_fp, index_col = 0)\n",
    "        compare_results = pd.concat([lme4_results, sm_results, mp_results], join=\"outer\", axis = 1)\n",
    "        compare_results[\"lme4/mp % diff\"] = \\\n",
    "                (compare_results[\"lme4\"] - compare_results[\"MPlus\"])/compare_results[\"MPlus\"]      \n",
    "        compare_results[\"sm/mp % diff\"] = \\\n",
    "                (compare_results[\"statsmodels\"] - compare_results[\"MPlus\"])/compare_results[\"MPlus\"]\n",
    "\n",
    "    # If we only want to compare R and Python models\n",
    "    else:\n",
    "        compare_results = pd.concat([lme4_results, sm_results], join=\"outer\", axis = 1)\n",
    "        compare_results[\"sm/lme4 % diff\"] = \\\n",
    "                (compare_results[\"lme4\"] - compare_results[\"statsmodels\"])/compare_results[\"statsmodels\"]\n",
    "\n",
    "    return compare_results\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ## How do lme4 (R) and Statsmodels (Python) compare when using REML?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will load test blood-pressure data and fit a simple mixed effects models in `R` and `Python`.`SEQN` is the grouping variable (subject's unique ID). [Based on this example.](https://github.com/kshedden/Statsmodels-MixedLM/blob/master/nhanes_bp.py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify `lme4` model implicitly using `REML`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "library(lme4)\n",
    "library(tidyverse)\n",
    "library(MplusAutomation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "bp_data = read.csv(\"bp_data.csv\")\n",
    "r_reml_model = lmer(bp ~ BMXBMI + female + (1 | SEQN), data = bp_data)\n",
    "lme4_params(r_reml_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will replicate this lme4 model with statsmodels. Again, we will be using `REML` estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify and fit model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "bp_data = pd.read_csv(\"bp_data.csv\")\n",
    "sm_reml_model = sm.MixedLM.from_formula(\"bp ~ female + BMXBMI\", groups = \"SEQN\", data = bp_data)\n",
    "sm_reml_results = sm_reml_model.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lme4</th>\n",
       "      <th>statsmodels</th>\n",
       "      <th>sm/lme4 % diff</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Int</th>\n",
       "      <td>97.535339</td>\n",
       "      <td>97.535</td>\n",
       "      <td>3.47952e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMXBMI</th>\n",
       "      <td>0.799370</td>\n",
       "      <td>0.799</td>\n",
       "      <td>0.000463643</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female</th>\n",
       "      <td>-3.296478</td>\n",
       "      <td>-3.296</td>\n",
       "      <td>0.00014512</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Int_SE</th>\n",
       "      <td>2.989652</td>\n",
       "      <td>2.99</td>\n",
       "      <td>-0.000116445</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMXBMI_SE</th>\n",
       "      <td>0.106120</td>\n",
       "      <td>0.106</td>\n",
       "      <td>0.00113113</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female_SE</th>\n",
       "      <td>1.546270</td>\n",
       "      <td>1.546</td>\n",
       "      <td>0.000174376</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEQN_Var</th>\n",
       "      <td>270.933125</td>\n",
       "      <td>270.918</td>\n",
       "      <td>5.58297e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Resid_Var</th>\n",
       "      <td>15.805451</td>\n",
       "      <td>15.8059</td>\n",
       "      <td>-2.74691e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEQN_SD</th>\n",
       "      <td>16.460046</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Resid_SD</th>\n",
       "      <td>3.975607</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEQN_SE</th>\n",
       "      <td>NaN</td>\n",
       "      <td>5.582</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 lme4 statsmodels sm/lme4 % diff\n",
       "Int         97.535339      97.535    3.47952e-06\n",
       "BMXBMI       0.799370       0.799    0.000463643\n",
       "female      -3.296478      -3.296     0.00014512\n",
       "Int_SE       2.989652        2.99   -0.000116445\n",
       "BMXBMI_SE    0.106120       0.106     0.00113113\n",
       "female_SE    1.546270       1.546    0.000174376\n",
       "SEQN_Var   270.933125     270.918    5.58297e-05\n",
       "Resid_Var   15.805451     15.8059   -2.74691e-05\n",
       "SEQN_SD     16.460046         NaN            NaN\n",
       "Resid_SD     3.975607         NaN            NaN\n",
       "SEQN_SE           NaN       5.582            NaN"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reml_comparison = compare_params(\"r_df.csv\", sm_reml_results)\n",
    "reml_comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference between the two models is relatively small. The diffferences that do arise are, for the most part, due to the inspecific rounding of `statsmodels` and the increased specificity of `lme4`. We cannot include `MPlus` in this analysis of `REML` models because the software doesn't have a `REML` estimator. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `NaN` values result in each software package using differing techniques to calculate variation among grouping variable effects. This will occur again in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How do all three software compare when using ML (Maximum Likelihood) estimation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's see how well `statsmodels`,`lme4` and `MPlus` line up if they all use `ML` estimation. We will follow the same process as above except this time, we will fit an `MPlus` model as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "r_ml_model = lmer(bp ~ BMXBMI + female + (1 | SEQN), REML = FALSE, data = bp_data)\n",
    "lme4_params(r_ml_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "sm_ml_model = sm.MixedLM.from_formula(\"bp ~ female + BMXBMI\", groups = \"SEQN\", data = bp_data)\n",
    "sm_ml_results = sm_ml_model.fit(reml = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write `MPlus` `ML` model script to hard drive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "echo \"\n",
    "TITLE:\n",
    "    Replication of Nhanes bp models;\n",
    "\n",
    "DATA:\n",
    "    FILE = \n",
    "        mplus_bp_data.csv;\n",
    "\n",
    "VARIABLE:\n",
    "    NAMES = !might need to be ARE instead of !\n",
    "        bp\n",
    "        RIDAGEYR\n",
    "        female\n",
    "        BMXBMI\n",
    "        SEQN;\n",
    "\n",
    "    USEVARIABLES = \n",
    "        bp\n",
    "        BMXBMI\n",
    "        female\n",
    "        SEQN;\n",
    "\n",
    "    CLUSTER = \n",
    "        SEQN;\n",
    "\n",
    "    WITHIN = \n",
    "        BMXBMI;\n",
    "        \n",
    "    BETWEEN = \n",
    "        female;\n",
    "\n",
    "ANALYSIS:\n",
    "    TYPE = TWOLEVEL RANDOM;\n",
    "    ESTIMATOR = ML\n",
    "\n",
    "MODEL:\n",
    "    %BETWEEN%\n",
    "        bp ON \n",
    "        female;\n",
    "    %WITHIN%\n",
    "        bp ON \n",
    "        BMXBMI;\n",
    "\" >| bp_ml.inp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses the [Mplus free demo version](https://www.statmodel.com/demo.shtml). There are options below for running this on Linux and MacOS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Linux\n",
    "#export PATH=$PATH:/opt/mplusdemo  \n",
    "#export TMP=/tmp\n",
    "#mpdemo bp_ml.inp   >/dev/null\n",
    "\n",
    "# MacOS\n",
    "/Applications/MplusDemo/mpdemo bp_ml.inp  >/dev/null"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract parameters from MPlus output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Reading model:  bp_ml.out \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R\n",
    "\n",
    "mplus_params(\"bp_ml.out\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lme4</th>\n",
       "      <th>statsmodels</th>\n",
       "      <th>MPlus</th>\n",
       "      <th>lme4/mp % diff</th>\n",
       "      <th>sm/mp % diff</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Int</th>\n",
       "      <td>97.535308</td>\n",
       "      <td>97.535</td>\n",
       "      <td>97.530</td>\n",
       "      <td>0.000054</td>\n",
       "      <td>5.12663e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMXBMI</th>\n",
       "      <td>0.799371</td>\n",
       "      <td>0.799</td>\n",
       "      <td>0.799</td>\n",
       "      <td>0.000464</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female</th>\n",
       "      <td>-3.296465</td>\n",
       "      <td>-3.296</td>\n",
       "      <td>-3.297</td>\n",
       "      <td>-0.000162</td>\n",
       "      <td>-0.000303306</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Int_SE</th>\n",
       "      <td>2.980047</td>\n",
       "      <td>2.98</td>\n",
       "      <td>2.980</td>\n",
       "      <td>0.000016</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BMXBMI_SE</th>\n",
       "      <td>0.105779</td>\n",
       "      <td>0.106</td>\n",
       "      <td>0.106</td>\n",
       "      <td>-0.002085</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>female_SE</th>\n",
       "      <td>1.541302</td>\n",
       "      <td>1.541</td>\n",
       "      <td>1.541</td>\n",
       "      <td>0.000196</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEQN_Var</th>\n",
       "      <td>269.161014</td>\n",
       "      <td>269.147</td>\n",
       "      <td>269.194</td>\n",
       "      <td>-0.000123</td>\n",
       "      <td>-0.000174595</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Resid_Var</th>\n",
       "      <td>15.805492</td>\n",
       "      <td>15.8059</td>\n",
       "      <td>15.806</td>\n",
       "      <td>-0.000032</td>\n",
       "      <td>-6.31062e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEQN_SD</th>\n",
       "      <td>16.406127</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Resid_SD</th>\n",
       "      <td>3.975612</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEQN_SE</th>\n",
       "      <td>NaN</td>\n",
       "      <td>5.535</td>\n",
       "      <td>17.955</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-0.691729</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Resid_SE</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.733</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 lme4 statsmodels    MPlus  lme4/mp % diff sm/mp % diff\n",
       "Int         97.535308      97.535   97.530        0.000054  5.12663e-05\n",
       "BMXBMI       0.799371       0.799    0.799        0.000464            0\n",
       "female      -3.296465      -3.296   -3.297       -0.000162 -0.000303306\n",
       "Int_SE       2.980047        2.98    2.980        0.000016            0\n",
       "BMXBMI_SE    0.105779       0.106    0.106       -0.002085            0\n",
       "female_SE    1.541302       1.541    1.541        0.000196            0\n",
       "SEQN_Var   269.161014     269.147  269.194       -0.000123 -0.000174595\n",
       "Resid_Var   15.805492     15.8059   15.806       -0.000032 -6.31062e-06\n",
       "SEQN_SD     16.406127         NaN      NaN             NaN          NaN\n",
       "Resid_SD     3.975612         NaN      NaN             NaN          NaN\n",
       "SEQN_SE           NaN       5.535   17.955             NaN    -0.691729\n",
       "Resid_SE          NaN         NaN    0.733             NaN          NaN"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ml_comparison = compare_params(\"r_df.csv\", sm_ml_results, \"mplus_df.csv\")\n",
    "ml_comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In an Mplus twolevel random model, there are 4 `ESTIMATOR` options: `MLR`, `ML`, `MLF`, and `Bayes`. `MLF` and `Bayes` produce standard errors much different than `R` and `Python`. `ML` `MPlus` estimation produces results that are very close to `Python` and `R` so that is where the overlap occurs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There isn't anything to compare the MPlus `MLR` model because neither `Python` nor `R` have an `MLR` option."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In R, there is a package called `robustlmm` which has a `rlmer` model that features two more estimators. Neither of these estimators are as close as `lme4::lmer`, though."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <span style=\"text-decoration: underline\">Conclusion</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using Maximum Likelihood estimation, all three software produce relatively congruent models. However, `R` and `Python` don't have any other methods that produce models similar to those of `MPlus`"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}