标签: mixed-models

我目前正在编写一个脚本来评估用于线性混合模型的(受限制的)对数似然函数.我需要它来计算模型的可能性,其中一些参数固定为任意值.也许这个脚本对你们中的一些人也有帮助!

我用lmer()从lme4和logLik()来检查我的脚本是否正常工作,因为它应该.而且看起来,它没有!由于我的教育背景并不真正关注这一级别的统计数据,我有点迷失了.

接下来,您将找到一个使用sleepstudy-data的简短示例脚本:

  # * * * * * * * * * * * * * * * * * * * * * * * * 
  # * example data

  library(lme4)
  data(sleepstudy)
  dat <- sleepstudy[ (sleepstudy$Days %in% 0:4) & (sleepstudy$Subject %in% 331:333) ,]
  colnames(dat) <- c("y", "x", "group")

  mod0 <- lmer( y ~ 1 + x + ( 1 | group ), data = dat)  


  # + + + + + + …

我发现了mlogit- 多项logit模型的包,用于寻找多项式混合logit模型.在阅读了优秀的插图后,我发现我无法在任何描述的示例中应用我的数据.

我现在写信希望对我的问题有所帮助,并创建了一个简单的例子来说明我的情况.

问题如下:在某处有辅音'Q'的单词.现在,我们进行了一项实验,他们的任务是听取这些话,并说他们是否听过Q,U或其他辅音.这必须依赖于诸如音节位置或真实/非真实单词之类的一些因素来建模.

在最小的例子中,我用音节位置创建了4个人和他们的答案.

library(mlogit)
library(nnet)
set.seed(1234)
data <- data.frame(personID = as.factor(sample(1:4, 40, replace=TRUE)),
               decision = as.factor(sample(c("Q","U", "other"), 40, replace=TRUE)),
               syllable = as.factor(sample(1:4, 40, replace=TRUE)))
summary(data)
 personID  decision  syllable
 1:11     other:10   1:18    
 2:10     Q    :18   2: 9    
 3:10     U    :12   3: 5    
 4: 9                4: 8 

据我所知nnet,multinom功能不包括混合型号.

modNnet1 <- multinom(decision ~ syllable, data=data)

首先,我使用mlogit.data-function来重塑文件.在与同事讨论后,我们得出结论,没有替代品.特定.变量.

 dataMod <- mlogit.data(data, shape="wide", choice="decision", id.var="personID")

 mod1 <- mlogit(formula = decision ~ 0|syllable,
           data = …

我想lme4::lmer在以下格式的60M观测值上运行混合模型(使用); 所有预测变量/因变量都是连续因变量的分类(因子)tc; patient是随机拦截项的分组变量.我有64位R和16Gb RAM,我在中央服务器上工作.RStudio是最新的服务器版本.

model <- lmer(tc~sex+age+lho+atc+(1|patient),
              data=master,REML=TRUE)

lho sex tc      age         atc patient
18  M   16.61   45-54       H   628143
7   F   10.52   12-15       G   2013855
30  M   92.73   35-44       N   2657693
19  M   24.92   70-74       G   2420965
12  F   17.44   65-69       A   2833610
31  F   7.03    75 and over A   1090322
3   F   28.59   70-74       A   2718649
29  F   4.09    75 and over C   384578
16  F   67.22   65-69       R   1579355
23  F   7.7 …

弗兰克·哈雷尔（Frank Harrell）的R包rms是实现多重逻辑回归的绝佳工具。但是，我想知道如何/是否有可能将随机效应纳入通过均方根的模型中。我知道rms可以通过nlme运行，但是只能通过广义最小二乘函数（Gls）运行，而不能通过lme函数运行，因此可以纳入随机效应。混合效应模型对于分析/解释可能会出现问题，但有时是必需的，以便考虑模型中的嵌套效应。

我不确定在这种情况下是否有用，但是我已经从rms帮助文件中复制了一些代码，这些文件运行了简单的逻辑回归模型，并添加了一行代码，表明通过MASS软件包的glmmPQL运行的混合效应逻辑回归模型。

n <- 1000    # define sample size
require(rms)
set.seed(17) # so can reproduce the results
age            <- rnorm(n, 50, 10)
blood.pressure <- rnorm(n, 120, 15)
cholesterol    <- rnorm(n, 200, 25)
sex            <- factor(sample(c('female','male'), n,TRUE))
label(age)            <- 'Age'      # label is in Hmisc
label(cholesterol)    <- 'Total Cholesterol'
label(blood.pressure) <- 'Systolic Blood Pressure'
label(sex)            <- 'Sex'
units(cholesterol)    <- 'mg/dl'   # uses units.default in Hmisc
units(blood.pressure) <- 'mmHg'
ch <- cut2(cholesterol, g=40, levels.mean=TRUE) # use mean values in intervals
table(ch) …

我有以下混合效应模型的输出。我想谈谈模型解释了多少变化。随机效应下的方差是否对应于残差（注：这里的试验是随机效应）所解释的变异？即 58.6 % 或者有其他方法来推断这一点

REML criterion at convergence: 71.9\n\nScaled residuals: \n     Min       1Q   Median       3Q      Max \n-1.82579 -0.59620  0.04897  0.62629  1.54639 \n\nRandom effects:\n Groups   Name        Variance Std.Dev.\n trial     (Intercept) 0.06008  0.2451  \n Residual              0.58633  0.7974  \nNumber of obs: 60, groups:  trial, 30\n\nFixed effects:\n                        Estimate Std. Error      df t value Pr(>|t|)    \n(Intercept)               1.5522     0.2684 12.6610  13.233  0.09888 \ndrugantho                 0.8871     0.1753 14.0000   1.043  0.31601    \ninterventionadded         0.2513     0.2553 14.0000  -1.276  0.32436 **   \nsexmale                   3.0026     0.6466 15.0000   4.066  0.00021  \n---\nSignif. codes:  0 \xe2\x80\x98***\xe2\x80\x99 …

我正在尝试获取混合模型预测的置信区间。预测函数不输出任何置信区间。很少有 StackOverflow 答案建议使用 merTools 包中的 PredictInterval 函数来获取间隔，但是这两个函数的预测估计之间存在差异，我试图在下图中进行比较。有人可以让我知道我在这里做错了什么吗？另外，我尝试构建的实际模型与下面代码片段中显示的模型类似，其中除了截距之外我没有固定效果组件。

library(merTools)
library(lme4)

dat <- iris
mod <- lmer(Sepal.Length ~ 1 + (1 + Sepal.Width + Petal.Length + 
                                  Petal.Width|Species), data=dat)

c1 <- predict(mod, dat)
c2 <- predictInterval(mod, dat)

plot_data <- cbind(c1, c2)
plot_data$order <- c(1:nrow(plot_data))

library(ggplot2)
ggplot(plot_data) + geom_line(aes(x=order, y=c1), color='red') + 
  geom_ribbon(aes(x=order, ymin=lwr, ymax=upr), color='blue', alpha=0.2) +  
  geom_line(aes(x=order, y=fit), color='blue')

红线表示预测“c1”，蓝线表示预测“c2”

我试图理解 Python statsmodel 包提供的混合线性模型的结果。我想避免数据分析和解释中的陷阱。问题在数据加载/输出代码块之后。

加载数据并拟合模型：

import statsmodels.api as sm
import statsmodels.formula.api as smf
data = sm.datasets.get_rdataset("dietox", "geepack").data
md = smf.mixedlm("Weight ~ Time", data, groups=data["Pig"])
mdf = md.fit()
print mdf.summary()

Mixed Linear Model Regression Results
========================================================
Model:            MixedLM Dependent Variable: Weight    
No. Observations: 861     Method:             REML      
No. Groups:       72      Scale:              11.3669   
Min. group size:  11      Likelihood:         -2404.7753
Max. group size:  12      Converged:          Yes       
Mean group size:  12.0                                  
--------------------------------------------------------
             Coef.  Std.Err.    z    P>|z| [0.025 0.975]
--------------------------------------------------------
Intercept    15.724    0.788  19.952 0.000 14.179 …

我有跨国面板数据，我想知道IV对二元学生水平结果的影响DV

我想包括一个嵌套的随机效应，该效应考虑到学生所在的学校会影响结果，以及不同国家的学校存在显着差异：(1|country/school). 所以我开始的模型是：

model = glmer(DV ~ IV + (1|country/school), data=data, family = 'binomial')

我还想考虑时间趋势。起初我认为我应该做年份固定效应，但这些国家的政治发展随着时间的推移而有很大差异，我想抓住这一点，虽然 1991 年可能让 A 国的学校陷入混乱，但 1991 年可能是教育资金的好年头在 B 国。因此我认为我可能应该包括一个国家年固定效应，如下所示：

model = glmer(DV ~ IV + (1|country/school) + as.factor(country_year),
                data=data, family = 'binomial')

模型的随机效应是：

Random effects:
 Groups          Name        Variance  Std.Dev. 
 school:country (Intercept) 5.703e-02 2.388e-01
 country         (Intercept) 4.118e-15 6.417e-08
Number of obs: 627, groups:  school:country, 51; country, 22

当模型中已经包含国家随机效应时，包含国家年固定效应是否不正确？

问这个问题的另一种方法：我应该如何应对可能与事实school的一个子集country，并且country_year是的一个子集country，但也school还是country_year彼此的子集？

我正在尝试为stan_glmer.nb( rstanarm) 输出创建一个表，但是model_parameters从包中parameters抛出一个奇怪的错误，我不确定如何解决。也许这是一个错误。

sessionInfo()版本信息的缩短输出：

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

parameters_0.8.2
rstanarm_2.21.1

一个可重现的例子：

library(rstanarm)
library(parameters)

x<-rnorm(500)
dat<-data.frame(x=x,z=rep(c("A","B","C","D","E"),100), y=.2+x*.7)

mod1<-stan_glmer(y~x+(x|z),data=dat)

model_parameters(mod1, effects="all")

我将在此处为您节省输出，因为它并不重要，但该功能有效。现在负二项式模型：

dat.nb<-data.frame(x=rnorm(500),z=rep(c("A","B","C","D","E"),100),
                y=rnbinom(500,size=1,prob = .5))

mod2<-stan_glmer.nb(y~x+(x|z),data=dat.nb)

model_parameters(mod2, effects="all")

现在出现错误消息：

Error in `$<-.data.frame`(`*tmp*`, "parameter", value = c("(Intercept)",  : 
  replacement has 3 rows, data has 1

尽管在parameters0.10.1 版本中，@BenBolker 得到一个空白输出，而不是错误（见评论）。不管怎样，这个函数似乎不适用于rstanarm离散分布（见评论）。据我在帮助文档中看到，没有任何内容表明需要指定负二项式模型。此外，该功能适用​​于lme4模型：

library(lme4)
mod1<-lmer(y~x+(x|z),data=dat)

model_parameters(mod1, effects="all")

mod2<-glmer.nb(y~x+(x|z),data=dat.nb)

model_parameters(mod2, effects="all") …

例子

library(glmmTMB)
library(ggeffects)

## Zero-inflated negative binomial model
(m <- glmmTMB(count ~ spp + mined + (1|site),
              ziformula=~spp + mined, 
              family=nbinom2, 
              data=Salamanders, 
              na.action = "na.fail"))
summary(m)

ggemmeans(m, terms="spp")
spp   | Predicted |       95% CI
--------------------------------
GP    |      1.11 | [0.66, 1.86]
PR    |      0.42 | [0.11, 1.59]
DM    |      1.32 | [0.81, 2.13]
EC-A  |      0.75 | [0.37, 1.53]
EC-L  |      1.81 | [1.09, 3.00]
DES-L |      2.00 | [1.25, 3.21]
DF    |      0.99 | [0.61, 1.62]

ggeffects::ggeffect(m, terms="spp")
spp …