"""R lme4 wrapper"""
import numbers
import warnings
import numpy as np
from .. import definitions as dfn
from ..rtdc_dataset.core import RTDCBase
from .rlibs import rpy2
from . import rsetup
[docs]class Lme4InstallWarning(UserWarning):
pass
[docs]class Rlme4(object):
def __init__(self, model="lmer", feature="deform"):
"""Perform an R-lme4 analysis with RT-DC data
Parameters
----------
model: str
One of:
- "lmer": linear mixed model using lme4's ``lmer``
- "glmer+loglink": generalized linear mixed model using
lme4's ``glmer`` with an additional a log-link function
via the ``family=Gamma(link='log'))`` keyword.
feature: str
Dclab feature for which to compute the model
"""
#: modeling method to use (e.g. "lmer")
self.model = None
#: dclab feature for which to perform the analysis
self.feature = None
#: list of [RTDCBase, column, repetition, chip_region]
self.data = []
#: model function
self.r_func_model = "feature ~ group + (1 + group | repetition)"
#: null model function
self.r_func_nullmodel = "feature ~ (1 + group | repetition)"
self.set_options(model=model, feature=feature)
# Make sure that lme4 is available
if not rsetup.has_lme4():
warnings.warn("Installing lme4, this may take a while!",
Lme4InstallWarning)
rsetup.install_lme4()
rsetup.import_lme4()
[docs] def add_dataset(self, ds, group, repetition):
"""Add a dataset to the analysis list
Parameters
----------
ds: RTDCBase
Dataset
group: str
The group the measurement belongs to ("control" or
"treatment")
repetition: int
Repetition of the measurement
Notes
-----
- For each repetition, there must be a "treatment" and a
"control" ``group``.
- If you would like to perform a differential feature analysis,
then you need to pass at least a reservoir and a channel
dataset (with same parameters for `group` and `repetition`).
"""
assert group in ["treatment", "control"]
assert isinstance(ds, RTDCBase)
assert isinstance(repetition, numbers.Integral)
region = ds.config["setup"]["chip region"]
# make sure there are no doublets
for ii, dd in enumerate(self.data):
if dd[1] == group and dd[2] == repetition and dd[3] == region:
raise ValueError("A dataset with group '{}', ".format(group)
+ "repetition '{}', and ".format(repetition)
+ "'{}' region has already ".format(region)
+ "been added (index {})!".format(ii))
self.data.append([ds, group, repetition, region])
[docs] def check_data(self):
"""Perform sanity checks on ``self.data``"""
# Check that we have enough data
if len(self.data) < 3:
msg = "Linear mixed effects models require repeated " \
+ "measurements. Please add more repetitions."
raise ValueError(msg)
[docs] def fit(self, model=None, feature=None):
"""Perform (generalized) linear mixed-effects model fit
The response variable is modeled using two linear mixed effect
models:
- model :const:`Rlme4.r_func_model` (random intercept +
random slope model)
- the null model :const:`Rlme4.r_func_nullmodel` (without
the fixed effect introduced by the "treatment" group).
Both models are compared in R using "anova" (from the
R-package "stats" :cite:`Everitt1992`) which performs a
likelihood ratio test to obtain the p-Value for the
significance of the fixed effect (treatment).
If the input datasets contain data from the "reservoir"
region, then the analysis is performed for the differential
feature.
Parameters
----------
model: str (optional)
One of:
- "lmer": linear mixed model using lme4's ``lmer``
- "glmer+loglink": generalized linear mixed model using
lme4's ``glmer`` with an additional log-link function
via ``family=Gamma(link='log'))`` :cite:`lme4`
feature: str (optional)
dclab feature for which to compute the model
Returns
-------
results: dict
Dictionary with the results of the fitting process:
- "anova p-value": Anova likelyhood ratio test (significance)
- "feature": name of the feature used for the analysis
``self.feature``
- "fixed effects intercept": Mean of ``self.feature`` for all
controls; In the case of the "glmer+loglink" model, the intercept
is already backtransformed from log space.
- "fixed effects treatment": The fixed effect size between the mean
of the controls and the mean of the treatments relative to
"fixed effects intercept"; In the case of the "glmer+loglink"
model, the fixed effect is already backtransformed from log
space.
- "fixed effects repetitions": The effects (intercept and
treatment) for each repetition. The first axis defines
intercept/treatment; the second axis enumerates the repetitions;
thus the shape is (2, number of repetitions) and
``np.mean(results["fixed effects repetitions"], axis=1)`` is
equivalent to the tuple (``results["fixed effects intercept"]``,
``results["fixed effects treatment"]``) for the "lmer" model.
This does not hold for the "glmer+loglink" model, because
of the non-linear inverse transform back from log space.
- "is differential": Boolean indicating whether or not
the analysis was performed for the differential (bootstrapped
and subtracted reservoir from channel data) feature
- "model": model name used for the analysis ``self.model``
- "model converged": boolean indicating whether the model
converged
- "r anova": Anova model (exposed from R)
- "r model summary": Summary of the model (exposed from R)
- "r model coefficients": Model coefficient table (exposed from R)
- "r stderr": errors and warnings from R
- "r stdout": standard output from R
"""
self.set_options(model=model, feature=feature)
self.check_data()
# Assemble dataset
if self.is_differential():
# bootstrap and compute differential features using reservoir
features, groups, repetitions = self.get_differential_dataset()
else:
# regular feature analysis
features = []
groups = []
repetitions = []
for dd in self.data:
features.append(self.get_feature_data(dd[1], dd[2]))
groups.append(dd[1])
repetitions.append(dd[2])
# Fire up R
with rsetup.AutoRConsole() as ac:
r = rpy2.robjects.r
# Load lme4
rpy2.robjects.packages.importr("lme4")
# Concatenate huge arrays for R
r_features = rpy2.robjects.FloatVector(np.concatenate(features))
_groups = []
_repets = []
for ii in range(len(features)):
_groups.append(np.repeat(groups[ii], len(features[ii])))
_repets.append(np.repeat(repetitions[ii], len(features[ii])))
r_groups = rpy2.robjects.StrVector(np.concatenate(_groups))
r_repetitions = rpy2.robjects.IntVector(np.concatenate(_repets))
# Register groups and repetitions
rpy2.robjects.globalenv["feature"] = r_features
rpy2.robjects.globalenv["group"] = r_groups
rpy2.robjects.globalenv["repetition"] = r_repetitions
# Create a dataframe which contains all the data
r_data = r["data.frame"](r_features, r_groups, r_repetitions)
# Random intercept and random slope model
if self.model == 'glmer+loglink':
r_model = r["glmer"](self.r_func_model, r_data,
family=r["Gamma"](link='log'))
r_nullmodel = r["glmer"](self.r_func_nullmodel, r_data,
family=r["Gamma"](link='log'))
else: # lmer
r_model = r["lmer"](self.r_func_model, r_data)
r_nullmodel = r["lmer"](self.r_func_nullmodel, r_data)
# Anova analysis (increase verbosity by making models global)
# Using anova is a very conservative way of determining
# p values.
rpy2.robjects.globalenv["Model"] = r_model
rpy2.robjects.globalenv["NullModel"] = r_nullmodel
r_anova = r("Anova = anova(Model, NullModel)")
try:
pvalue = r_anova.rx2["Pr(>Chisq)"][1]
except ValueError: # rpy2 2.9.4
pvalue = r_anova[7][1]
r_model_summary = r["summary"](r_model)
r_model_coefficients = r["coef"](r_model)
try:
fe_reps = np.array(r_model_coefficients.rx2["repetition"])
except ValueError: # rpy2 2.9.4
fe_reps = np.concatenate((
np.array(r_model_coefficients[0][0]).reshape(1, -1),
np.array(r_model_coefficients[0][1]).reshape(1, -1)),
axis=0)
r_effects = r["data.frame"](r["coef"](r_model_summary))
try:
fe_icept = r_effects.rx2["Estimate"][0]
fe_treat = r_effects.rx2["Estimate"][1]
except ValueError: # rpy2 2.9.4
fe_icept = r_effects[0][0]
fe_treat = r_effects[0][1]
if self.model == "glmer+loglink":
# transform back from log
fe_treat = np.exp(fe_icept + fe_treat) - np.exp(fe_icept)
fe_icept = np.exp(fe_icept)
fe_reps[:, 1] = np.exp(fe_reps[:, 0] + fe_reps[:, 1]) \
- np.exp(fe_reps[:, 0])
fe_reps[:, 0] = np.exp(fe_reps[:, 0])
# convergence
try:
lme4l = r_model_summary.rx2["optinfo"].rx2["conv"].rx2["lme4"]
except ValueError: # rpy2 2.9.4
lme4l = r_model_summary[17][3][1]
if lme4l and "code" in lme4l.names:
try:
conv_code = lme4l.rx2["code"]
except ValueError: # rpy2 2.9.4
conv_code = lme4l[0]
else:
conv_code = 0
ret_dict = {
"anova p-value": pvalue,
"feature": self.feature,
"fixed effects intercept": fe_icept,
"fixed effects treatment": fe_treat, # aka "fixed effect"
"fixed effects repetitions": fe_reps,
"is differential": self.is_differential(),
"model": self.model,
"model converged": conv_code == 0,
"r anova": r_anova,
"r model summary": r_model_summary,
"r model coefficients": r_model_coefficients,
"r stderr": ac.get_warnerrors(),
"r stdout": ac.get_prints(),
}
return ret_dict
[docs] def get_differential_dataset(self):
"""Return the differential dataset for channel/reservoir data
The most famous use case is differential deformation. The idea
is that you cannot tell what the difference in deformation
from channel to reservoir is, because you never measure the
same object in the reservoir and the channel. You usually just
have two distributions. Comparing distributions is possible
via bootstrapping. And then, instead of running the lme4
analysis with the channel deformation data, it is run with
the differential deformation (subtraction of the bootstrapped
deformation distributions for channel and reservoir).
"""
features = []
groups = []
repetitions = []
# compute differential features
for grp in sorted(set([dd[1] for dd in self.data])):
# repetitions per groups
grp_rep = sorted(set([dd[2] for dd in self.data if dd[1] == grp]))
for rep in grp_rep:
feat_cha = self.get_feature_data(grp, rep, region="channel")
feat_res = self.get_feature_data(grp, rep, region="reservoir")
bs_cha, bs_res = bootstrapped_median_distributions(feat_cha,
feat_res)
# differential feature
features.append(bs_cha - bs_res)
groups.append(grp)
repetitions.append(rep)
return features, groups, repetitions
[docs] def get_feature_data(self, group, repetition, region="channel"):
"""Return array containing feature data
Parameters
----------
group: str
Measurement group ("control" or "treatment")
repetition: int
Measurement repetition
region: str
Either "channel" or "reservoir"
Returns
-------
fdata: 1d ndarray
Feature data (Nans and Infs removed)
"""
assert group in ["control", "treatment"]
assert isinstance(repetition, numbers.Integral)
assert region in ["reservoir", "channel"]
for dd in self.data:
if dd[1] == group and dd[2] == repetition and dd[3] == region:
ds = dd[0]
break
else:
raise ValueError("Dataset for group '{}', repetition".format(group)
+ " '{}', and region".format(repetition)
+ " '{}' not found!".format(region))
fdata = ds[self.feature][ds.filter.all]
fdata_valid = fdata[~np.logical_or(np.isnan(fdata), np.isinf(fdata))]
return fdata_valid
[docs] def is_differential(self):
"""Return True if the differential feature is computed for analysis
This effectively just checks the regions of the datasets
and returns True if any one of the regions is "reservoir".
See Also
--------
get_differential_features: for an explanation
"""
for dd in self.data:
if dd[3] == "reservoir":
return True
else:
return False
[docs] def set_options(self, model=None, feature=None):
"""Set analysis options"""
if model is not None:
assert model in ["lmer", "glmer+loglink"]
self.model = model
if feature is not None:
assert dfn.scalar_feature_exists(feature)
self.feature = feature