标签: mixed-models

我用geepackR来估算逻辑边际模型geeglm().但我得到垃圾估计.它们大约16个数量级太大.然而,p值似乎与我的预期相似.这意味着响应基本上成为阶梯函数.见附图

以下是生成图表的代码:

require(geepack)
data = read.csv(url("http://folk.uio.no/mariujon/data.csv"))
fit = geeglm(moden ~ 1 + power, id = defacto, data=data, corstr = "exchangeable", family=binomial)
summary(fit)
plot(moden ~ power, data=data)
x = 0:2500
y = predict(fit, newdata=data.frame(power = x), type="response" )
lines(x,y)

这是回归表:

Call:
geeglm(formula = moden ~ 1 + power, family = binomial, data = data, 
    id = defacto, corstr = "exchangeable")

 Coefficients:
             Estimate   Std.err  Wald Pr(>|W|)    
(Intercept) -7.38e+15  1.47e+15  25.1  5.4e-07 ***
power        2.05e+13  1.60e+12 164.4  < 2e-16 …

我有一些重复的措施，顺序响应数据：

dat <- data.frame(
  id = factor(sample(letters[1:5], 50, replace = T)),
  response = factor(sample(1:7, 50, replace = T), ordered = T),
  x1 = runif(n = 50, min = 1, max = 10),
  x2 = runif(n = 50, min = 100, max = 1000)
)

我建立了以下模型：

library(ordinal)

model <- clmm(response ~ x1 + x2 + (1|id), data = dat)

我有一些新数据：

new_dat <- data.frame(
  id = factor(sample(letters[1:5], 5, replace = T)),
  x1 = runif(n = 5, min = 1, max = 10), …

假设我们假设了以下线性混合模型（LMM），我们通常将其称为fit。

library(lme4)

fit <- lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy)

假设我们也有兴趣从 中提取原始公式（即固定效应和随机效应部分）fit，我们希望将其传递给名为 的第二个 LMM fit2。显然，传递terms(fit)给formula的参数lmer()是行不通的。

> fit2 <- lmer(terms(fit), data = sleepstudy)
Error: No random effects terms specified in formula

问题

有没有办法从 中提取固定效应和随机效应部分fit，然后将其传递给lmer()？

关于更复杂的实验设计的混合模型有几个问题和帖子,所以我认为这个更简单的模型将有助于其他初学者以及我.

所以,我的问题是我想从sas proc混合程序中在R中制定重复测量ancova:

proc mixed data=df1;
FitStatistics=akaike
class GROUP person day;
model Y = GROUP X1 / solution alpha=.1 cl;
repeated / type=cs subject=person group=GROUP;
lsmeans GROUP;
run;

以下是使用R(下面)中创建的数据的SAS输出:

.           Effect       panel    Estimate       Error      DF    t Value    Pr > |t|     Alpha       Lower       Upper
            Intercept              -9.8693      251.04       7      -0.04      0.9697       0.1     -485.49      465.75
            panel        1         -247.17      112.86       7      -2.19      0.0647       0.1     -460.99    -33.3510
            panel        2               0           .       .        .         .             .           .           .
            X1                     20.4125     10.0228       7       2.04      0.0811       0.1      1.4235 …

我正在使用 lme4 包进行线性混合效果建模

混合效应模型如下：

fm01 <- lmer(sublat <- goal + (1|userid))

上面的命令返回一个名为 fm01 的 S4 对象

此对象包括固定效应及其 OLS 标准误差（如下）

Fixed effects:

            Estimate Std. Error t value
(Intercept)   31.644      3.320   9.530
goaltypeF1    -4.075      3.243  -1.257
goaltypeF2    -9.187      5.609  -1.638
goaltypeF3   -13.935      9.455  -1.474
goaltypeF4   -20.219      8.196  -2.467
goaltypeF5   -12.134      8.797  -1.379"

但是，我需要提供可靠的标准错误

如何使用 lme4 返回的 S4 对象执行此操作？

我正在使用MCMCglmm包中的一些贝叶斯线性混合模型R.我的数据包括使用错误测量的预测变量.因此,我想建立一个考虑到这一点的模型.我的理解是,基本的混合效应模型MCMCglmm将仅对响应变量(如ols回归)中的误差最小化.换句话说,垂直误差将被最小化.我想最小化与回归线/平面/超平面正交的误差.

1. 是否可以使用MCMCglmm或者我必须使用JAGS/STAN来进行变量误差(也就是总最小二乘)模型？
2. 是否可以在同一模型中使用多个预测变量来实现这一点(我有一些模型有3个或4个预测变量,每个模型都有误差测量)？
3. 如果可能,我该如何指定模型？

我在下面列出了一个数据集,其中一个随机变量height用错误来衡量,以说明基本设置MCMCglmm.

library(nlme)
library(MCMCglmm)

data(Orthodont)

set.seed(1234)

Orthodont$height <- c(rnorm(54, 170, 10), rnorm(54, 150, 10))

prior1 <- list(
    B = list(mu = rep(0, 3), V = diag(1e+08, 3)), 
    G = list(G1 = list(V = 1, nu = 1, alpha.mu = 0, alpha.V = 1000)), 
    R = list(V = 1, nu = 0.002)
)

model1 <- MCMCglmm(
    fixed = distance ~ height + Sex,   
    random …

我发现该predict函数目前尚未在使用R包clmm中的函数拟合的累积链接混合模型中实现ordinal。虽然在同一个包中predict实现clmm2，但我选择应用clmm，因为后者允许多个随机效果。此外，我还安装了几个clmm模型并使用包model.avg中的函数执行模型平均MuMIn。理想情况下，我想使用平均模型来预测概率。然而，虽然MuMIn支持clmm模型，predict也不适用于普通模型。

有没有办法破解该predict函数，以便该函数不仅可以预测clmm模型的概率，还可以使用模型平均系数clmm（即“平均”类的对象）进行预测？例如：

require(ordinal)
require(MuMIn)

mm1 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
        link = "probit", threshold = "equidistant")

## test random effect:
mm2 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|RESP:PROD), data = soup,
        link = "logistic", threshold = "equidistant")

#create a model selection object
mm.sel<-model.sel(mm1,mm2) …

我正在尝试将该glm方法caret::train用于广义线性混合效应模型。我使用以下代码来设置该功能：

`GLMERmod <- list(type="Classification", library="lme4", loop=NULL)

parameters <- data.frame(parameter="parameter", class="character", 
label="parameter")

GLMERmod$parameters <- parameters


grid <- function (x, y, len = NULL, search = "grid"){
  data.frame(parameter = "none")
 } 



 GLMERmod$grid <- grid


 fit <- 
 function (x, y, wts, param, lev, last, classProbs, ...) 
 {
    dat <- if (is.data.frame(x)) 
        x
    else as.data.frame(x)
    dat$.outcome <- y
    if (length(levels(y)) > 2) 
        stop("glm models can only use 2-class outcomes")
    theDots <- list(...)
    if (!any(names(theDots) == "family")) {
         theDots$family <- if …