标签: mixed-models

我有一个混合模型，数据如下所示：

> head(pce.ddply)
  subject Condition errorType     errors
1    j202         G         O 0.00000000
2    j202         G         P 0.00000000
3    j203         G         O 0.08333333
4    j203         G         P 0.00000000
5    j205         G         O 0.16666667
6    j205         G         P 0.00000000

每个受试者为 errorType（O 或 P）提供两个数据点，并且每个受试者处于条件 G（N=30）或 N（N=33）。errorType 是一个重复变量，Condition 是一个变量之间。我对主效应和相互作用都感兴趣。所以，首先是一个方差分析：

> summary(aov(errors ~ Condition * errorType + Error(subject/(errorType)),
                 data = pce.ddply))

Error: subject
          Df  Sum Sq  Mean Sq F value Pr(>F)
Condition  1 0.00507 0.005065   2.465  0.122
Residuals 61 0.12534 0.002055               

Error: subject:errorType
                    Df  Sum …

我已经构建了一个lme4用于模型选择的模型，dredge但是我无法将随机效应与相关的固定效应对齐。我的完整模型的结构如下。

fullModel<-glmer(y ~x1 + x2 + (0+x1|Year) + (0+x1|Country) + (0+x2|Year) + (0+x2|Country) + (1 | Year) +(1|Country), family=binomial('logit'),data = alldata)

在这个模型结构中，模型选择dredge产生了三种固定效应的组合，即 x1、x2 和 x1+x2，但是随机效应结构保持与完整模型相同，这样即使固定效应只有 x1，随机效应将包括(0+x2|Year) + (0+x2|Country). 例如，只有 x1 作为固定效应的模型，在随机效应结构中仍然会有 x2，如下所示。

y ~x1 + (0+x1|Year) + (0+x1|Country) + (0+x2|Year) +(0+x2|Country) + (1 | Year) +(1|Country), family=binomial('logit')

有没有办法配置dredge不选择其中指定了其他固定效果的随机效果？我有大约 x1….x50。

正如标题所示，我正在尝试从 lme4 （或其他包？）对象中提取矩阵。为了明确我想要什么，我认为最简单的方法是参考 SAS 文档：https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect022.htm

随机效应的方差-协方差矩阵

在 SAS 表示法中，该矩阵称为G ，是随机效应参数gamma的方差-协方差矩阵。通过使用 PROC MIXED 和输出传输系统中的选项“G”，您可以获得G作为矩阵。

我知道一旦我有了伽玛的方差分量和维度，手动构建这个矩阵就相对简单。尽管如此，我还是希望有一种更简单的方法。

混合模型方程解

在 SAS 表示法中，这些称为C。通过使用 PROC MIXED 和输出传送系统中的选项“MMEQSOL”，您可以请求生成混合模型方程的解以及逆系数矩阵。我感兴趣的是后者。

提前致谢！

import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv('http://www.bodowinter.com/tutorial/politeness_data.csv')
df = df.drop(38)

在R我会做：

lmer(frequency ~ attitude + (1|subject) + (1|scenario), data=df)

这R给了我：

Random effects:
 Groups   Name        Variance Std.Dev.
 scenario (Intercept)  219     14.80   
 subject  (Intercept) 4015     63.36   
 Residual              646     25.42   
Fixed effects:
            Estimate Std. Error t value
(Intercept)  202.588     26.754   7.572
attitudepol  -19.695      5.585  -3.527

我试图用做同样的事情statsmodels：

model = smf.mixedlm("frequency ~ attitude", data=df, groups=df[["subject","scenario"]]).fit()

但是model.summary()给了我不同的输出：

      Mixed Linear Model Regression Results
=======================================================
Model: …

我对连续结果进行了重复测量，如下所示：

library(magrittr)
library(ggplot2)
library(nlme)
mydata <- structure(list(ID = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 
7L, 7L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 11L, 11L, 
12L, 12L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 16L, 16L, 
17L, 17L, 17L, 17L, 17L, 18L, 18L, 19L, 19L, 20L, 20L, 21L, 21L, 
22L, 22L, 22L, 22L, 22L, 23L, 23L, 24L, …

我想获得模型交互效果的预测（例如使用sjPlot::plot_model或ggeffects打包，然后将它们提供给 ggplot2 以可视化 ggplot2 中的交互。有人可以帮助编写代码来做到这一点吗？我的问题与堆栈上的其他问题不同，因为它首先使用标准包进行预测，然后绘制它们。

我当前的型号：

model <- glmmTMB(total_count ~ mean_temp*lwd_duration + (1|year), family = nbinom1, data=df)\n\nsummary(model)\n\nFamily: nbinom1  ( log )\nFormula:          total_count ~ mean_temp * lwd_duration + (1 | year)\nData: df\n\n     AIC      BIC   logLik deviance df.resid \n   260.6    270.7   -124.3    248.6       34 \n\nRandom effects:\n\nConditional model:\n Groups Name        Variance  Std.Dev. \n year   (Intercept) 3.683e-09 6.069e-05\nNumber of obs: 40, groups:  year, 4\n\nDispersion parameter for nbinom1 family ():  178 \n\nConditional model:\n                        Estimate Std. Error z value Pr(>|z|) …

当我尝试使用 glm 运行我的模型时出现此错误。为什么 R 不使用我的嵌套字符嵌套效果？

model <- glm(Propotion ~ Parasites_box * Parasites_nest + Day_of_year + Site + (1|Subsite/Nestbox/Cavity_ID), data = data.ex, family = binomial)

Error in Subsite/Trapnest : non-numeric argument to binary operator

使用 glmer 给了我一个不同的警告：

model <- glmer(Propotion ~ Parasites_box * Parasites_nest + Day_of_year + Site + (1|Subsite/Nestbox/Cavity_ID), data = data.ex, family = binomial)

Error in (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L, maxit = 100L,  : 
  Downdated VtV is not positive definite
In addition: …