Strutural Equation Modelling (SEM) with Jupyter and Anaconda¶
This notebook is intended as an example of how to "do SEM" using different tools that are accessible in a Jupyter IPython Notebook installed with Anaconda. The notebook uses a Python3 kernel and interfaces with several R packages. Alternatives would be to run an R kernel in the notebook directly, or to use an RStudio markdown notebook and the reticulate package to access Python functionality from R.
To Do for future revisions:
Integrate instructions on using different hosted / collaborative notebook options for this tutorial:
- CoCalc: https://cocalc.com/
- Google Colaboratory in Google Drive: https://colab.research.google.com/
- (Python works, conda & R aren't available so most of the meat of this notebook won't work. If you block third-party cookies, you might have to allow cookies from [*.]googleusercontent.com to get output working)
- Microsoft Azure Notebooks: https://notebooks.azure.com/
- Binder: https://mybinder.org/
- A JupyterHub server: https://jupyterhub.readthedocs.io/en/stable/ (if I ever get to installing one)
Turn it into a RISE presentation https://rise.readthedocs.io/en/docs_hot_fixes/ (ideally also using the "Split Cells Notebook" extension for better visual organization)
Integrate other SEM tools available in R:
- lavaan, openmx, semplot (DONE)
- sem (less well documented/supported)
- metasem and metafor for meta-analysis (goes beyond normal SEM)
- simsem's semtools
- bigsem, ctsem https://github.com/cdriveraus/ctsem, gsem, lavaan.survey, regsem, semdiag
Compare with other tools:
- Onyx GUI for SEM http://onyx.brandmaier.de/
For local use, I recommend jupyterlab as the interface, usually acessible as http://localhost:8888/lab
Installing Required (and Optional) Components¶
To run the code in this notebook, a local (or hosted) installation of Anaconda is required. The following installs R components available through Anaconda. We use the --yes option to avoid prompts for confirmation (impossible in the notebook) but that also means you cannot check what will be installed before it proceeds. To avoid this, you could run the installations from an Anaconda prompt or using the Anaconda Navigator.
R Studio is not strictly required, but it is a useful alternative way of using R. Also note that this assumes an Anaconda installation/environment that you have permission to change. On Windows, that means choosing the recommended "for this user only" option on install of Anaconda. Otherwise, these installations require Administrator / root rights.
!conda install --yes r-essentials
# !conda install --yes rstudio
!conda install --yes r-lavaan
Installing and Enabling the Python-R bridge¶
Now we need rpy2 which provides the bridge between this python notebook and R. (Also installing tzlocal due to a current dependency bug in rpy2.)
!conda install --yes tzlocal rpy2
Let's enable the rpy2 extension, so that we can then execute R code with the %%R magic command at the top of a cell.
%load_ext rpy2.ipython
If you are on Windows, and you do not get text output from %%R cells showing up in the notebook, but instead in the console window where jupyter is running, this is a bug in rpy2 on Windows, and there's a workaround to capture stdout by running the following cells, see https://github.com/vitorcurtis/RWinOut
%%R
install.packages(c("R.utils"))
!curl -O "https://raw.githubusercontent.com/vitorcurtis/RWinOut/master/RWinOut.py"
%load_ext RWinOut
Installing Interface Conveniences¶
And finally we want to see plots directly in the notebook. The simplest way is to request inline plots.
%matplotlib inline
But assuming you are running this in a JuyterLab interface, you might want the ipympl library to get interactive widget plots.
# We install ipympl with pip, as it is not yet readily available with conda:
!pip install ipympl
# nodejs is needed for the interactive features if using JupyterLab,
# the corresponding package for normal notebooks, widgetsnbextension, should already be installed.
!conda install --yes nodejs
# install the jupyterlab extensions:
!jupyter labextension install @jupyter-widgets/jupyterlab-manager
!jupyter labextension install jupyter-matplotlib
Now enable widget-based plots:
%matplotlib widget
Finally, since this is a long document, the following extension adds a table-of-contents sidebar to the JupyterLab interface:
!jupyter labextension install @jupyterlab/toc
Installing More R Packages¶
Now for some R packages that are not readily available in an Anaconda default install. They might be available through the conda-forge "channel" - however, at the time of writing, I cannot recommend this, as the performance of conda install is abysmal when using R packages from that repository.
When using a server-hosted notebook, some or all of these packages might already be installed.
%%R
install.packages(c("semPlot", "OpenMx", "semTools", "sem", "gpairs", "GGally"))
A note to avoid possible confusion: lavaan provides a function cfa as a convenience for confirmatory factor analysis. There is also an R package called cfa - however, that one is not related to SEM.
Importing and Configuring Python & R Packages¶
First, import all the basic packages in Python, such as pandas, numpy, matplotlib.
Also import seaborn for simple high-level plots with decent looks. It complements the default plotting provided by matplotlib. See here, for a useful brief overview of looking at data using some of seaborn's plot types: https://elitedatascience.com/python-seaborn-tutorial
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
Next, some of the libraries that are statistics-oriented.
import statsmodels.api as sm
import statsmodels.formula.api as smf # the R-like interface for statsmodels
import statsmodels.graphics as smg
import sklearn
Now for loading R packages, we will be using.
%%R
library(lavaan)
library(semPlot)
library(OpenMx)
library(semTools)
In R, there are also several packages providing convenient high-level plots, such as generalized pairs plots.
%%R
library(GGally)
library(gpairs)
Configure Plotting¶
For the purpose of getting nice visual output, we will also set some defaults for plotting libraries.
defaultfigwidth, defaultfigheight = 10, 9
# set a slightly larger default size for plots. default dpi is 100
plt.rcParams['figure.figsize'] = [defaultfigwidth, defaultfigheight]
# enable seaborn's defaults for nicer plots overall:
sns.set(color_codes=True)
There is currently no built-in way to set default dimensions for plots in R. The %R magic command from the rpy2 library accepts width, height, and units parameters like this: %%R -w 10 -h 9 -u in -r 100 but it would be nice to set defaults.
Since this is python, there is a way around that using monkey-patching. Note that this is usually a Bad Idea(TM) and should be avoided if possible. It is also purely cosmetic for the purposes of this notebook, so it can be safely ignored. :)
# these are the defaults we want to set:
default_units = 'in' # inch, to make it more easily comparable to matpplotlib
default_res = 100 # dpi, same as default in matplotlib
default_width = 10
default_height = 9
# try monkey-patching a function in rpy2, so we effectively get these
# default settings for the width, height, and units arguments of the %R magic command
import rpy2
old_setup_graphics = rpy2.ipython.rmagic.RMagics.setup_graphics
def new_setup_graphics(self, args):
if getattr(args, 'units') is not None:
if args.units != default_units: # a different units argument was passed, do not apply defaults
return old_setup_graphics(self, args)
args.units = default_units
if getattr(args, 'res') is None:
args.res = default_res
if getattr(args, 'width') is None:
args.width = default_width
if getattr(args, 'height') is None:
args.height = default_height
return old_setup_graphics(self, args)
rpy2.ipython.rmagic.RMagics.setup_graphics = new_setup_graphics
Loading Data¶
We are borrowing example data from this excellent course offered at Harvard, S090A1: https://canvas.harvard.edu/courses/8737/pages/data
The actual data is from the Zambian Early Childhood Development Project. The full sample has more than 1600 Zambian six-year-olds, from a study led by Günther Fink and Stephanie Zuilkowski.
The data is in the proprietory stata format, so we first need to convert it and import it. We will use pandas.read_stata but this could also be accomplished in R with the foreign package. First, we write a helper function to download the data file if it is not in the current directly. Defining a function will allow us to re-use it later for other datasets.
import os.path
import urllib.request
def downloadIfMissing(filenameData, remoteLocation):
'''Check if the file exists. If not, try downloading from remoteLocation.'''
if not os.path.isfile(filenameData):
with urllib.request.urlopen(remoteLocation) as response:
with open(filenameData, 'xb') as destinationFile:
destinationFile.write(response.read())
# make sure we have the small data file available in the current directory, if not, try to download it:
filenameSmallZambiaData = "S090_InClass_Zambia.dta"
downloadIfMissing(filenameSmallZambiaData, "https://canvas.harvard.edu/courses/8737/files/1839865/download")
# read the data into a pandas dataframe
smallZambiaDF = pd.read_stata(filenameSmallZambiaData)
len(smallZambiaDF) # should return 1613
# make sure we have the full measurement data file available in the current directory, if not, try to download it:
filenameMeasureZambiaData = "S090_InClass_Zambia_Measurement.dta"
downloadIfMissing(filenameMeasureZambiaData, "https://canvas.harvard.edu/courses/8737/files/1994882/download")
# read the data into a pandas dataframe
measureZambiaDF = pd.read_stata(filenameMeasureZambiaData)
len(measureZambiaDF) # should return 1623
Data Overview¶
We now have the data in a dataframe, let's get an overview of what kind of data we are dealing with.
Now let's have a look at the dataframe.
smallZambiaDF.head() # equivalent to [:5] i.e. first five entries
smallZambiaDF.info()
smallZambiaDF.describe(include='all', percentiles=[]) # describe categorical and numerical columns, don't bother with percentiles
smallZambiaDF.male.describe()
smallZambiaDF.male.value_counts()
smallZambiaDF.wealth.value_counts()
# visually check relations between numeric variables
grid = sns.pairplot(smallZambiaDF, hue="ece", height=defaultfigheight/6, kind='scatter')
# we have to use an explicit height per facet-figure here, since a grid of figures doesn't follow the matplotlib default size
grid = sns.pairplot(smallZambiaDF, hue="ece", height=defaultfigheight/6, kind='reg') # linear regressions on top of scatter
pd.scatter_matrix(smallZambiaDF)
plt.show()
plt.figure()
sns.swarmplot(x="male", y="socemo", data=smallZambiaDF)
plt.show()
%%R -i smallZambiaDF
pairs(smallZambiaDF)
%%R -i smallZambiaDF
ggpairs(smallZambiaDF) # from library GGally