标签: mixed-models

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

所以,我的问题是我想从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 …

在R的nlme包中的lme()函数的标准示例中:

fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
summary(fm2)

出现了相关表:

Correlation: 
          (Intr) age   
age       -0.813       
SexFemale -0.372  0.000

如果涉及许多因素组合,这可能是巨大的.

有没有办法在summary命令中抑制输出？我知道我可以用

   print(fm2, cor=F) 

但这并没有向我显示通常输出的其余部分,例如没有p值计算.

我在R中对森林的特定处理进行了荟萃分析.对于这个模型,我需要拟合随机效应来解释方法的研究差异和站点年龄的变化之间,因为这两者都是混淆变量,我对调查由它们引起的变化没有明确的兴趣.

但是,据我所知,[metfor]当你有一个多级模型时,包不允许你计算一个R平方类型的统计量.

无论如何,在这里更清楚地描述我的问题是一个模拟数据集

Log<-data.frame(Method=rep(c("RIL","Conv"),each=10),
     RU=runif(n=20,min=10,max=50),SDU=runif(n=20,5,20),
     NU=round(runif(n=20,10,20),0))
Log$Study<-rep(1:4,each=5)
Log$Age<-rep(c(0,10,15,10),times=5)
RIL<-(Log$RU-(Log$RU*(abs(rnorm(n=20,mean=.6,sd=0.1)))))+(0.5*Log$Age)
Conv<-(Log$RU-(Log$RU*(abs(rnorm(n=20,mean=.2,sd=0.1)))))+(0.2*Log$Age)
Log$RL<-ifelse(Log$Method=="RIL",RIL,Conv)
Log$SDL<-Log$SDU
Log$NL<-Log$NU

#now we perform a meta-analysis using metafor
require(metafor)
ROM<-escalc(data=Log,measure="ROM",m2i=RU,
sd2i=SDU,n2i=NU,m1i=RL,sd1i=SDL,n1i=NL,append=T)
Model1<-rma.mv(yi,vi,random=~(1|Study)+(1|Age),method="ML",data=ROM)
summary(Model1)
forest(Model1)

上述模型是一个空模型,用于观察截距在统计上是否与零显着不同.就我们而言.但是,我还想看看治疗方面的差异是否描述了我在森林情节中看到的效果大小的差异,你可以在这里看到

所以我运行这个模型:

Model2<-rma.mv(yi,vi,mods=~Method,random=~(1|Study)+(1|Age),method="ML",data=ROM)
summary(Model2)

哪个看起来不错.

    Multivariate Meta-Analysis Model (k = 20; method: ML)

  logLik  Deviance       AIC       BIC      AICc  
  0.4725   19.8422    7.0550   11.0380    9.7217  

Variance Components: 

outer factor: Age   (nlvls = 3)
inner factor: Study (nlvls = 4)

            estim    sqrt  fixed
tau^2      0.0184  0.1357     no
rho        1.0000             no

Test for Residual Heterogeneity: 
QE(df …

我正在使用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 …

我想在模型中使用nlme::lme(底部的数据)指定不同的随机效果.随机效应是:1)intercept并且position变化subject; 2) intercept变化comparison.这很简单lme4::lmer:

lmer(rating ~ 1 + position + 
     (1 + position | subject) + 
     (1 | comparison), data=d)

> ...
Random effects:
 Groups     Name        Std.Dev. Corr 
 comparison (Intercept) 0.31877       
 subject    (Intercept) 0.63289       
            position    0.06254  -1.00
 Residual               0.91458      
 ...

但是,我想坚持,lme因为我也想模拟自相关结构(position是一个时间变量).我怎么能像上面那样使用lme？我在下面的尝试会影响效果,这不是我想要的.

lme(rating ~ 1 + position,
random = list( ~ 1 + position | subject,
               ~ 1 | comparison), data=d) …

我有这种格式的数据集(df)

index <- runif(n = 100,min = 0, max = 1)
type1 <- rep("low", 50)
type2 <- rep("high", 50)
type <- c(type1,type2)

level1 <- rep("single", 25)
level2 <- rep("multiple", 25)
level3 <- rep("single", 25)
level4 <- rep("multiple", 25)
level <- c(level1,level2,level3,level4)

block <- rep(1:5, 10)
set <- rep(1:5, 10)

df <- data.frame("index" = index,"type" = type, "level" = level, "block" = block, "set" = set)
df$block <- as.factor(df$block)
df$set <- as.factor(df$set)

我想创建一个看起来像这样的模型

model <- lmer(index ~ type * level + …

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

我用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)  


  # + + + + + + …

传统上，线性混合效应模型是按以下方式制定的。Ri = Xi×？+ Zi×bi +？i在哪里？代表估计的固定效应，Z代表随机效应。因此，X是经典设计矩阵。使用R，在使用nlme包中的lme拟合模型后，我希望能够提取这两个矩阵。例如，也可以在nlme软件包中找到的数据集“ Rails”包含在6条随机选择的铁路轨道上的三个超声行进时间的单独测量值。我可以通过以下方法为每个轨道设置截距固定效果和随机效果的简单模型。

library(nlme)
lmemodel<-lme(travel ~ 1, random = ~ 1 | Rail, data=Rail)

X设计矩阵只是一个18x1矩阵（6个滑轨* 3个测量值）的单位，并且可以通过以下方式轻松提取：

 model.matrix(lmemodel, data=Rail)
   (Intercept)
1            1
2            1
3            1
4            1
5            1
6            1
7            1
8            1
9            1
10           1
11           1
12           1
13           1
14           1
15           1
16           1
17           1
18           1
attr(,"assign")
[1] 0

我想做的是提取随机效果设计矩阵Z。我意识到，如果我使用lme4包拟合相同的模型，则可以通过以下方式完成：

library(lme4)
lmermodel<-lmer(travel ~ 1 + (1|Rail),data=Rail) 
t(lmermodel@Zt)  ##takes the transpose of lmermodel@Zt
lmermodel@X  ## extracts the X matrix …

我有以下混合效应模型的输出。我想谈谈模型解释了多少变化。随机效应下的方差是否对应于残差（注：这里的试验是随机效应）所解释的变异？即 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 …

我试图理解 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 …