标签: mixed-models

aov可以通过以下model.tables函数获得最小二乘平均值及其对象的标准误差:

npk.aov <- aov(yield ~ block + N*P*K, npk)
model.tables(npk.aov, "means", se = TRUE)

我想知道如何得到广义最小二乘意味着它们的标准误差nlme或lme4对象:

library(nlme)
data(Machines)
fm1Machine <- lme(score ~ Machine, data = Machines, random = ~ 1 | Worker )

任何评论和提示将受到高度赞赏.谢谢

4个蜂箱配备了传感器,可收集蜂巢内的温度,湿度,压力和分贝.这些是响应变量.

治疗是wifi暴露,实验组从第1天到第20天暴露于wifi,然后从第35-45天再次暴露,并且收集数据直到第54天.荨麻疹数= 4,n由传感器收集的数据每个蜂巢=〜百万.

我在运行混合效果模型时遇到困难.

所有蜂房的响应变量都有一个数据框.

names(Hives)
[1] "time"           "dht22_t"        "dht11_t"        "dht22_h"       
[5] "dht11_h"        "db"             "pa"             "treatment_hive"
[9] "wifi"   

时间为"%Y-%m-%d%H:%M:%S",dht11/22_t/h是温度和湿度数据."wifi"是对应于暴露时间的二分变量(1 =开0 =关),治疗蜂巢是暴露于wifi的荨麻疹的另一个二分变量(1 =暴露,0 =对照).

这是我得到的错误.

attach(Hives)
model2 = lme(pa_t~wifi*treatment_hive, random=time, na.action=na.omit, method="REML",)

Error in reStruct(random, REML = REML, data = NULL) : 
Object must be a list or a formula

以下是代码示例:

    time    dht22_t dht11_t dht22_h dht11_h db  pa  treatment_hive  wifi
1   01/09/2014 15:19    NA      NA  NA      NA  51.75467    NA      0   1
2   01/09/2014 15:19    30.8    31  59.8    44  55.27682    100672  0 …

在 R 中，我正在寻找一种方法来估计lmer使用 kenward-rogers 或 satterthwaite 自由度和 SE 的模型的线性对比的置信区间。

例如，我可以使用 t 值（来自 KR 的 df）和 SE 计算混合模型中固定效应参数的 CI，如带有 R 的 SAS。

mod<-lmerTest::lmer(y~time1+treatment+time1:treatment+(1|PersonID),data=data)
lmerTest::summary(mod,ddf = "Kenward-Roger")

这个输出：

Fixed effects:
                Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)      49.0768     1.0435 56.4700  47.029  < 2e-16 ***
time1             5.8224     0.5963 48.0000   9.764 5.51e-13 ***
treatment         1.6819     1.4758 56.4700   1.140   0.2592    
time1:treatment   2.0425     0.8433 48.0000   2.422   0.0193 * 

允许时间 1 的 CI，例如：

5.8224+abs(qt(0.05/2, 48))*0.5963 #7.021342
5.8224-abs(qt(0.05/2, 48))*0.5963 #4.623458

我想为固定系数的线性对比做同样的事情。这是 p 值，但没有 SE 输出。 …

是否有一种简单的方法来拟合 R 中的多元回归，其中因变量根据Skellam 分布（两个泊松分布计数之间的差异）分布？就像是：

myskellam <- glm(A ~ B + C + D, data = mydata, family = "skellam")

这应该适应固定效应。但理想情况下，我更喜欢随机效应，因为我知道固定效应可能会引入测量偏差。因此我想理想的解决方案应该是使用lme4orglmmADMB包。

或者，有没有办法转换数据以应用更常用的回归工具？

我知道这是一个有点老的问题，但我只是想知道现在是否有解决方案。

lme4我通常使用带有函数的包来执行混合效果模型lmer。但是，我知道这个函数不允许我在模型中包含负方差分量。

我真的很想在 R 模型中包含负方差。有没有人对我会使用哪些包有任何建议？或者说，新的lme4允许吗？

我目前正在阅读 Andy Field 的书《使用 R 发现统计》。第 14 章是关于混合建模的，他使用了包lme中的函数nlme。

他使用快速约会数据创建的模型是这样的：

speedDateModel <- lme(dateRating ~ looks + personality +
                 gender + looks:gender + personality:gender + 
                 looks:personality,
    random = ~1|participant/looks/personality)

我尝试使用包lmer中的函数重新创建一个类似的模型lme4；然而，我的结果不同。我以为我有正确的语法，但也许不是？

speedDateModel.2 <- lmer(dateRating ~ looks + personality + gender + 
              looks:gender + personality:gender + 
              (1|participant) + (1|looks) + (1|personality), 
              data = speedData, REML = FALSE)

另外，当我运行这些模型的系数时，我注意到它只为每个参与者产生随机截距。然后我试图创建一个可以产生随机截距和斜率的模型。我似乎无法为这两个函数提供正确的语法来执行此操作。任何帮助将不胜感激。

我使用包中的lmer()函数构建了一个混合效果模型lme4。由于lme4某些很好的哲学原因，该包不输出系数的 p 值。但是，我仍然需要在我的出版物中报告 p 值。我知道有多种方法可以使用由创建的模型来计算 p 值lmer()，例如这里。

我的问题是：我想使用包中的tidy()函数提取 p 值broom。在这里，我真的很想坚持，tidy()因为我想维护以下管道：

data_frame %>% group_by(grouping variables) %>% do(tidy(fitted_model))

一种选择是创建一个自定义函数并将其附加到管道中。但是，该broom软件包的手册页（http://rpackages.ianhowson.com/cran/broom/man/lme4_tidiers.html）说：

"p.value  P-value computed from t-statistic (may be missing/NA)". 

通过这个，我假设一个从 lmer 模型给出的 t 值计算 p 值的函数已经在 broom 中实现了。所以，我不愿意重新发明轮子。

问题是我根本没有得到名为 p.value 的列。我期待一个名为 p.value 的列，其中 NA 作为最坏的情况。

代码：

library(lme4)
library(broom)
lme <- lmer(Reaction ~ Days + (1 + Days | Subject), sleepstudy)
tidy(lme)
tidy(lme, effects = …

我正在尝试使用包中的lmer函数来拟合混合模型lme4。但是，我不明白应该向start参数输入什么。我的目的是使用简单的线性回归来使用估计的系数作为混合模型的起始值。

假设我的模型如下：

linear_model = lm(y ~ x1 + x2 + x3, data = data)
coef = summary(linear_model)$coefficients[- 1, 1] #I remove the intercept
result = lmer(y ~ x1 + x2 + x3 | x1 + x2 + x3, data = data, start = coef)

这个例子是我正在做的事情的一个过于简化的版本，因为我将无法共享我的数据。

然后我得到以下类型的错误：

Error during wrapup: incorrect number of theta components (!=105) #105 is the value I get from the real regression I am trying to fit.

我尝试了许多不同的解决方案，试图提供一个列表并命名这些值，theta就像我在某些论坛上看到的那样。

另外 …

我尝试使用R重现一些SAS输出.我想重现的方法是:

使用混合模型重复测量因子时间的方差的双向分析(协方差矩阵= CS,估计方法= REML)

AIC看起来很好看...我想知道是否有人知道SAS使用的AIC公式......

主要的SAS输出是:

如果loglik相同,Anova表是相同的,但不是AIC(和BIC)事件.

这就是我用R做的:

library(nlme)
dataset_melt <- structure(list(Groupe = c("A", "A", "A", "A", "A", "B", "B", 
"B", "B", "B", "C", "C", "C", "C", "C", "A", "A", "A", "A", "A", 
"B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "A", "A", "A", 
"A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", "C", "A", 
"A", "A", "A", "A", "B", "B", "B", "B", "B", "C", "C", "C", "C", 
"C", "A", "A", "A", "A", "A", "B", "B", …

输入数据,从https://pastebin.com/1f7VuBkx复制(太大,不包括在这里)

data.frame':    972 obs. of  7 variables:
$ data_mTBS : num  20.3 22.7 0 47.8 58.7 ...
$ data_tooth: num  1 1 1 1 1 1 1 1 1 1 ...
$ Adhesive  : Factor w/ 4 levels "C-SE2","C-UBq",..: 2 2 2 2 2 2 2 2 2 2 ...
$ Approach  : Factor w/ 2 levels "ER","SE": 1 1 1 1 1 1 1 1 1 1 ...
$ Aging     : Factor w/ 2 levels "1w","6m": 1 …