text: # coding: utf-8
  
  # DO NOT EDIT
  # Autogenerated from the notebook mixed_lm_example.ipynb.
  # Edit the notebook and then sync the output with this file.
  #
  # flake8: noqa
  # DO NOT EDIT
  
  # # 线性混合效应模型
  
  import numpy as np
  import statsmodels.api as sm
  import statsmodels.formula.api as smf
  
  # 比较 R 的 lmer 模型和 statsmodels 的 MixedLM 模型
  # =======================================
  #
  # statsmodels 的线性混合模型（MixedLM）的插值紧密遵循 Lindstrom 和 Bates（JASA 1988）中概述的方法。 
  # 这也是 R 软件包 LME4 遵循的方法。 其他软件包如Stata，SAS等，也遵循此方法，因为该领域的基本技术已很成熟。
  #
  # 在这里，我们将展示如何使用 statsmodels 的 MixedLM 拟合线性混合模型的过程，并与 R（LME4）的结果进行比较。
  # 这是我们的导入申明:
  
  # ## 猪的生长曲线
  #
  # 这些是析因实验的纵向数据。 结果变量是每头猪的体重，这里我们将使用的唯一预测变量是“时间”。 
  # 首先，我们拟合一个模型，该模型将平均体重表示为时间的线性函数，每头猪都有随机截距。
  # 使用公式指定模型。 由于未指定随机效应结构，因此会自动使用默认的随机效应结构（每个组的随机一个截距）。
  data = sm.datasets.get_rdataset('dietox', 'geepack').data
  md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"])
  mdf = md.fit()
  print(mdf.summary())
  
  # 这是使用 R 的 LMER 拟合的相同模型:
  
  # 请注意，在statsmodels 的 summary 结果中，固定效应和随机效应参数估计值显示在一个单表中。 
  # 在 statsmodels 上面输出的动物随机效应标记为 "Intercept RE"。 在 LME4 输出中，这个效应
  # 是在随机效应下猪的截距。
  #
  # 关于随机效应方差和协方差参数的标准误差是否有用存在很多争论。 在LME4中，不会显示这些标准误差，
  # 因为该软件包的作者认为它们不是非常有用。 尽管有充分的理由质疑其效用，但我们选择在 summary 表中
  # 包括标准误差，但未显示相应的 Wald 置信区间。
  #
  # 接下来，我们为每个动物拟合一个具有两种随机效应的模型：随机截距和随机斜率（相对于时间）。这意味着
  # 每头猪可能具有不同的基线体重，并且以不同的速度生长。该公式指定“时间”是具有随机系数的协变量。默认情况下，
  # 公式始终包含一个截距（此处可以使用 "0 + Time" 作为公式将其抑制）。
  
  md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"], re_formula="~Time")
  mdf = md.fit()
  print(mdf.summary())
  
  # 这是使用 R 的 LMER 拟合的相同模型::
  
  # $(0.294 / \sqrt{19.493 * 0.416} \approx 0.1)$ 的随机截距和随机斜率仅是弱相关的
  # 因此，接下来我们拟合一个模型，其中两个随机效应是不相关：
  
  
  .294 / (19.493 * .416)**.5
  
  md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"], re_formula="~Time")
  free = sm.regression.mixed_linear_model.MixedLMParams.from_components(
      np.ones(2), np.eye(2))
  
  mdf = md.fit(free=free)
  print(mdf.summary())
  
  # 当我们将相关性参数固定为 0 时，似然降低了 0.3。将 2 x 0.3 = 0.6 与 chi^2 1 df 参考分布进行比较表明，
  # 数据与参数等于 0 的模型非常一致。
  #
  # 
  # 这是使用 R 的 LMER 拟合的相同模型（请注意，此处 R 报告的是 REML 标准而不是似然，其中 REML标准 是对数似然的两倍）：Here is the same model fit using LMER in R (note that here R is
  
  
  # ## Sitka 生长数据
  #
  #这是 R 库的 LMER 提供的示例数据集之一。 结果变量是树的大小，此处使用的协变量是时间值。数据按树分组。
  
  data = sm.datasets.get_rdataset("Sitka", "MASS").data
  endog = data["size"]
  data["Intercept"] = 1
  exog = data[["Intercept", "Time"]]
  
  # 这是 statsmodels 的 LME 拟合的一个随机截距基本模型。 我们通过 Endog 和 Exog 直接作为数组传递给 LME init 函数。
  # 还要注意，endog_re 明确指定了参数4作为一个随机截距（尽管这可能是默认值，而不是指定）。
  
  md = sm.MixedLM(endog, exog, groups=data["tree"], exog_re=exog["Intercept"])
  mdf = md.fit()
  print(mdf.summary())
  
  # 这是使用 R 的 LMER 拟合的相同模型:
  
  # 现在我们可以尝试添加一个随机斜率。这次我们从R开始。从下面的代码和输出中，我们看到随机斜率方差的 REML 估计几乎为零。
  
  # 如果我们使用 statsmodels 的 LME 的默认设置运行此命令，则会看到方差估计的确很小，这会导致有关解决方案
  # 在参数空间的边界上的警告。回归斜率与 R (拟合模型的斜率)非常吻合，但是似然值远高于 R 返回的似然值。
  
  
  exog_re = exog.copy()
  md = sm.MixedLM(endog, exog, data["tree"], exog_re)
  mdf = md.fit()
  print(mdf.summary())
  
  # 我们可以通过绘制轮廓似然图来进一步探索随机效应结构。我们从随机截距开始，生成一个从 0.1个单位之下到0.1个单位之上的
  # MLE 的轮廓似然图。 因为配置文件内的每个优化都会产生警告（由于随机斜率方差接近于零），因此我们在此关闭警告。
  
  
  import warnings
  
  with warnings.catch_warnings():
      warnings.filterwarnings("ignore")
      likev = mdf.profile_re(0, 're', dist_low=0.1, dist_high=0.1)
  
  # 这是一个轮廓似然图函数。我们将对数似然差乘以2，使它服从带有 1个自由度的 $\chi^2$ 参考分布。
  
  
  import matplotlib.pyplot as plt
  
  plt.figure(figsize=(10, 8))
  plt.plot(likev[:, 0], 2 * likev[:, 1])
  plt.xlabel("Variance of random slope", size=17)
  plt.ylabel("-2 times profile log likelihood", size=17)
  
  # 这是一个轮廓似然图函数。轮廓似然图显示随机斜率方差参数的 MLE 模型是一个非常小的正数，并且此估计的不确定性较低。
  
  re = mdf.cov_re.iloc[1, 1]
  likev = mdf.profile_re(1, 're', dist_low=.5 * re, dist_high=0.8 * re)
  
  plt.figure(figsize=(10, 8))
  plt.plot(likev[:, 0], 2 * likev[:, 1])
  plt.xlabel("Variance of random slope", size=17)
  plt.ylabel("-2 times profile log likelihood", size=17)
  