标签: mixed-models

我正在使用lmer()包lme4来估计混合效果模型.这很有效,但现在我想在固定数量的迭代中运行估算过程,然后通过指定由上一个估算过程计算的起始值来恢复过程.

根据这方面的帮助,?lmer可以通过设置参数:

• start- 这些是新的起始值,根据帮助,可以ST从拟合模型中提取槽中的值并使用这些值,即使用x@ST
• maxiter - 作为命名参数提供给 control

因此,例如,假设我想要lme使用iris数据,可以尝试这样做:

library(lme4)

# Fit model with limited number of iterations

frm <- "Sepal.Length ~ Sepal.Width | Species"

x <- lmer(frm, data=iris, 
          verbose=TRUE, control=list(maxIter=1), model=FALSE)

# Capture starting values for next set of iterations
start <- list(ST=x@ST)

# Update model
twoStep <-  lmer(frm, data=iris, 
          verbose=TRUE, control=list(maxIter=100), model=TRUE, 
          start=start)

这有效.看一下输出,其中第一列是REML,即随机效应最大似然.特别注意模型2中的REML从模型1终止的地方开始:

> x <- lmer(frm, data=iris, 
+           verbose=TRUE, control=list(maxIter=1), …

我有来自参与者的受试者内部生理数据(part),他们都在三轮(round)中查看了刺激(阅读报纸),每轮都有五篇论文(paper),并且在报纸中每个都有不同的访问次数(visit) .我有两个固定的因素(CONDhier和CONDabund)加上相互作用来预测生理状态(例如EDA),这通常是自回归的.我试着考虑生理学中随机效应的个体差异(让我们暂时解决拦截问题),并且可能会因为其他随机效应而导致疲劳.

因此,我希望在R中运行的模型在SPSS中是:

MIXED EDA BY CONDhier CONDabund  
 /FIXED=CONDhier CONDabund CONDhier*CONDabund | SSTYPE(3)  
 /RANDOM=INTERCEPT | SUBJECT(part) COVTYPE(VC)  
 /RANDOM=INTERCEPT | SUBJECT(part*round) COVTYPE(VC)  
 /PRINT=SOLUTION  
 /METHOD=REML  
 /REPEATED=visit | SUBJECT(part*round*paper) COVTYPE(AR1).

现在,我已经明白,虽然lme没有做好交叉术语,lmer(处理交叉术语没有问题)不能使用不同的协方差结构.我可以运行简单的lme模型,如

    lme(EDA ~ factor(CONDhier) * factor(CONDabund), random= ~1
   |part, na.action=na.exclude, data=phys2)

但是更复杂的模型超出了我的范围.我已经读过lme中的交叉术语可以用随机定义来完成

    random=pdBlocked(list(pdCompSymm(~part), pdCompSymm(~round-1), pdCompSymm(~paper-1), 
pdCompSymm(~visit-1)))

但这似乎阻止了AR1结构,并且第二次随机拦截了部分*round,来自我.而且我不太确定它和我的SPSS语法一样.

那么,有什么建议吗？虽然在lme和lmer上有很多不同的着作,但我找不到一个既有交叉术语又有AR1的作品.

(另外,关于lme的语法看起来很模糊:从几个不同的来源我已经理解了|在左边的东西下面是什么,组成嵌套术语,~1是随机截距,~x是随机斜率,和~1 + x都是,但似乎至少有:和-1定义,我在任何地方都找不到.是否有一个教程可以解释所有不同的定义？)

在lm和glm模型中,我使用函数coef并confint实现目标:

m = lm(resp ~ 0 + var1 + var1:var2) # var1 categorical, var2 continuous
coef(m)
confint(m)

现在我添加了随机效果到模型 - 使用lmerlme4包的函数使用混合效果模型.但是,功能coef,confint不再为我工作!

> mix1 = lmer(resp ~ 0 + var1 + var1:var2 + (1|var3)) 
                                      # var1, var3 categorical, var2 continuous
> coef(mix1)
Error in coef(mix1) : unable to align random and fixed effects
> confint(mix1)
Error: $ operator not defined for this S4 class

我试图谷歌并使用文档,但没有结果.请指出我正确的方向.

编辑:我也在考虑这个问题是否更适合https://stats.stackexchange.com/但我认为它比统计更具技术性,所以我认为它最适合这里(SO)......你怎么看？

我有以下模型

x  <- rep(seq(0, 100, by=1), 10)
y  <- 15 + 2*rnorm(1010, 10, 4)*x + rnorm(1010, 20, 100)
id <- NULL
for(i in 1:10){ id <- c(id, rep(i,101)) }
dtfr <- data.frame(x=x,y=y, id=id)
library(nlme)
with(dtfr, summary(     lme(y~x, random=~1+x|id, na.action=na.omit)))
model.mx <- with(dtfr, (lme(y~x, random=~1+x|id, na.action=na.omit)))
pd       <- predict( model.mx, newdata=data.frame(x=0:100), level=0)
with(dtfr, plot(x, y))
lines(0:100, predict(model.mx, newdata=data.frame(x=0:100), level=0), col="darkred", lwd=7)

用predict,level=0我可以绘制平均人口反应.如何从nlme对象中提取和绘制整个群体的95%置信区间/预测带？

我可以使用nlme包中的gls()来构建没有随机效果的mod1.然后我可以将使用AIC的mod1与使用lme()构建的mod2进行比较,其中包含随机效果.

mod1 = gls(response ~ fixed1 + fixed2, method="REML", data)
mod2 = lme(response ~ fixed1 + fixed2, random = ~1 | random1, method="REML",data)
AIC(mod1,mod2)

是否有类似于gls()的lme4包,它允许我构建没有随机效果的mod3,并将它与使用lmer()构建的mod4进行比较,其中包含随机效果？

mod3 = ???(response ~ fixed1 + fixed2, REML=T, data)
mod4 = lmer(response ~ fixed1 + fixed2 + (1|random1), REML=T, data)
AIC(mod3,mod4)

我正在尝试运行混合效果模型,F2_difference该模型预测其余列作为预测变量,但我收到一条错误消息

固定效应模型矩阵排名不足,因此下降7列/系数.

从这个链接,固定效果模型是排名不足,我想我应该findLinearCombos在R包中使用caret.但是,当我尝试时findLinearCombos(data.df),它给了我错误信息

qr.default(object)中的错误:外部函数调用中的NA/NaN/Inf(arg 1)另外:警告消息:在qr.default(object)中:强制引入的NAs

我的数据没有任何NA - 可能导致这种情况的原因是什么？(对不起,如果答案很明显 - 我是R的新手).

我的所有数据都是除了我试图预测的数值之外的因素.这是我的数据的一小部分样本.

sex <- c("f", "m", "f", "m")
nasal <- c("TRUE", "TRUE", "FALSE", "FALSE")
vowelLabel <- c("a", "e", "i", "o")
speaker <- c("Jim", "John", "Ben", "Sally")
word_1 <- c("going", "back", "bag", "back")
type <- c("coronal", "coronal", "labial", "velar")
F2_difference <- c(345.6, -765.8, 800, 900.5)
data.df <- data.frame(sex, nasal, vowelLabel, speaker,
                      word_1, type, F2_difference
                      stringsAsFactors = TRUE)

编辑:这是一些更多的代码,如果它有帮助.

formula <- F2_difference ~ …

我需要在一个相对较大的数据框架上运行逻辑回归,其中包含480个条目和3个固定效果变量.固定效应var A有3233级,var B有2326级,var C有811级.总而言之,我有6370个固定效果.数据是横截面的.如果我不能使用正常glm函数运行此回归,因为回归矩阵对于我的记忆来说似乎太大了(我得到消息" Error: cannot allocate vector of size 22.9 Gb").我正在寻找在我的Macbook Air(OS X 10.9.5 8GB RAM)上运行此回归的替代方法.我也可以访问具有16GB RAM的服务器.

我试过用几种不同的方法解决这个问题,但到目前为止还没有取得令人满意的结果:

LFE/felm:使用的felm回归函数lfe包减去运行回归之前固定的效果.这完美地工作,并允许我在几分钟内将上述回归作为正常线性模型运行.但是,lfe不支持逻辑回归和glms.所以felm非常适合了解不同模型的模型拟合,但不适用于最终的逻辑回归模型.

biglm/bigglm:我想过bigglm用来将我的功能分解成更易于管理的块.然而,若干来源(例如link1,link2,link3)提到为了使其起作用,因子级别需要在块之间保持一致,即每个块必须包含每个因子变量的每个因子中的至少一个.因子A和B包含仅出现一次的级别,因此我无法将这些集合拆分为具有一致级别的不同块.如果我删除固定效应A的10个因子和B的8个因子(微小的变化),我将只剩下4个级别的因子,并且将我的数据分成4个块将使其更易于管理.然而,我仍然需要弄清楚如何对我的df进行排序,以确保我的480.000条目被分类为4个块,其中3个因子中的每个因子的每个因子级别出现至少一次.

GlmmGS/glmgs:glmmgs具有相同名称的包中的函数执行固定效果减法,如lfe使用"Gauss-Seidel"算法的逻辑回归包.不幸的是,该包已不再开发.对R来说比较新,没有深入的统计经验,我无法理解输出,也不知道如何以一种能给我正常"效果大小","模型拟合","模型拟合"的方式对其进行转换.显着性区间"glm回归摘要提供的指标.

我给包的作者发了一条消息.他们回应如下:

该包不提供与glm对象相同格式的输出.但是,在给定当前输出的情况下,您可以轻松计算大部分拟合统计量(估计的标准误差,拟合度)(在CRAN版本中,我相信当前输出是系数估计的向量,以及相关的向量标准误差;协方差分量相同,但如果你没有随机效应拟合模型,你不必担心它们).只要注意用于计算标准误差的协方差矩阵是与Gauss-Seidel算法相关的精度矩阵的对角线块的倒数,因此它们倾向于低估联合似然的标准误差.我不再维护包裹,我没有时间详细说明; 包装背后的开创性理论可以在手册中引用的论文中找到 ,其他一切都需要用笔和纸来制定:).

如果任何人都可以解释如何"轻松计算大部分拟合统计数据",使得没有任何统计学教育的人能够理解它(可能是不可能的)或者提供R代码,以示例如何实现这一点,我将是非常感谢!

Revolution Analytics:我在一台模拟Mac上的Windows 7的虚拟机上安装了革命分析企业.该程序具有一个被调用的函数RxLogit,该函数针对大型逻辑回归进行了优化.使用RxLogit我得到的功能the error (Failed to allocate 326554568 bytes. Error in rxCall("RxLogit", params) : bad allocation),所以该功能似乎也遇到了内存问题.但是,该软件使我能够在分布式计算集群上运行回归.所以我可以通过在具有大量内存的集群上购买计算时间来"解决问题".但是,我想知道革命分析程序是否提供了我不知道的任何公式或方法,这将允许我做某种类似的lfe固定效果减法操作或类似的bigglm …

在尝试R使用lme4包时改进线性混合效果模型的模型拟合时间时,我想分享一些我的想法.

数据集大小:数据集大约包含400,000行和32列.遗憾的是,无法分享有关数据性质的信息.

假设和检查:假设响应变量来自正态分布.在模型拟合过程之前,使用相关表和aliasR中提供的函数测试变量的共线性和多重共线性.

连续变量按比例缩放以帮助收敛.

模型结构:模型方程包含31个固定效应(包括截距)和30个随机效应(不包括截距).对于具有2700个级别的特定因子变量,随机效应是随机的.协方差结构是方差分量,因为假设随机效应之间存在独立性.

模型方程示例:

lmer(Response ~ 1 + Var1 + Var2 + ... + Var30 + (Var1-1| Group) + (Var2-1| Group) + ... + (Var30-1| Group), data=data, REML=TRUE)

模型已成功安装,但是,提供结果需要大约3.1小时.SAS中的相同型号花了几秒钟.网上有关于如何通过使用非线性优化算法减少时间的文献,以及在优化nloptwrap完成后执行的耗时衍生计算calc.derivs = FALSE:

https://cran.r-project.org/web/packages/lme4/vignettes/lmerperf.html

时间减少了78%.

问题:是否有其他替代方法可以通过相应地定义lmer参数输入来缩短模型拟合时间？在模型拟合时间方面,R和SAS之间存在很大差异.

任何建议表示赞赏.

我正在尝试在 python 中实现混合效应逻辑回归。作为比较，我使用的是 R 包中的glmer函数lme4。

我发现该statsmodels模块有一个BinomialBayesMixedGLM应该能够适合这样的模型。但是，我遇到了一些问题：

1. 我发现该statsmodels函数的文档并不完全有用或清晰，所以我不完全确定如何正确使用该函数。
2. 到目前为止，我的尝试还没有产生可以复制我glmer在 R 中拟合模型时得到的结果。
3. 我希望该BinomialBayesMixedGLM函数不计算 p 值，因为它是贝叶斯，但我似乎无法弄清楚如何访问参数的完整后验分布。

作为测试用例，我使用了可用的泰坦尼克号数据集here。

import os
import pandas as pd
import statsmodels.genmod.bayes_mixed_glm as smgb

titanic = pd.read_csv(os.path.join(os.getcwd(), 'titanic.csv'))

r = {"Pclass": '0 + Pclass'}
mod = smgb.BinomialBayesMixedGLM.from_formula('Survived ~ Age', r, titanic)
fit = mod.fit_map()
fit.summary()

#           Type    Post. Mean  Post. SD       SD SD (LB) SD (UB)
# Intercept M           3.1623    0.3616            
# Age       M …

我试图用纵向设计来估算数据集中的数据.有两个预测因子(实验组和时间)和一个结果变量(分数).聚类变量是id.

这是玩具数据

set.seed(345)
A0 <- rnorm(4,2,.5)
B0 <- rnorm(4,2+3,.5)
A1 <- rnorm(4,6,.5)
B1 <- rnorm(4,6+2,.5)
A2 <- rnorm(4,10,.5)
B2 <- rnorm(4,10+1,.5)
A3 <- rnorm(4,14,.5)
B3 <- rnorm(4,14+0,.5)
score <- c(A0,B0,A1,B1,A2,B2,A3,B3)
id <- rep(1:8,times = 4, length = 32)
time <- rep(0:3, each = 8, length = 32)
group <- rep(c("A","B"), times =2, each = 4, length = 32)
df <- data.frame(id = id, group = group, time = time,  score = score)

# plots
(ggplot(df, aes(x = time, y = …