标签: random-effects

qqmath函数使用lmer软件包的输出产生很好的随机效应图.也就是说,qqmath非常适合绘制层次模型中的截距及其在点估计周围的误差.下面使用名为Dyestuff的lme4包中的内置数据,给出了lmer和qqmath函数的一个示例.代码将使用ggmath函数生成分层模型和一个漂亮的图.

library("lme4")
data(package = "lme4")

# Dyestuff 
# a balanced one-way classiï¬cation of Yield 
# from samples produced from six Batches

summary(Dyestuff)             

# Batch is an example of a random effect
# Fit 1-way random effects linear model
fit1 <- lmer(Yield ~ 1 + (1|Batch), Dyestuff) 
summary(fit1)
coef(fit1) #intercept for each level in Batch 

# qqplot of the random effects with their variances
qqmath(ranef(fit1, postVar = TRUE), strip = FALSE)$Batch

最后一行代码产生了每个截距的非常好的图,每个估计周围都有误差.但格式化qqmath函数似乎非常困难,而且我一直在努力格式化情节.我想出了一些我无法回答的问题,我认为如果他们使用lmer/qqmath组合,其他人也可以从中受益:

1. 有没有办法采取上面的qqmath函数并添加一些选项,例如,使某些点为空而不是填充,或不同点的不同颜色？例如,您是否可以填充Batch变量的A,B和C的点数,但其余的点是否为空？
2. 是否可以为每个点添加轴标签(例如,可能沿顶部或右侧y轴)？
3. 我的数据接近45个截距,因此可以在标签之间添加间距,这样它们就不会相互碰撞？主要是,我感兴趣的是在图表上的点之间进行区分/标记,这在ggmath函数中似乎很麻烦/不可能.

到目前为止,在qqmath函数中添加任何附加选项会产生错误,如果它是标准图,我不会得到错误,所以我很茫然.

另外,如果你觉得有一个更好的包装/功能来绘制来自lmer输出的拦截,我很乐意听到它!(例如,你能用dotplot做点1-3吗？)

谢谢.

编辑:如果可以合理格式化,我也可以使用替代的dotplot.我只是喜欢ggmath情节的外观,所以我开始提出一个问题.

在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)......你怎么看？

我正在努力理解lmer的功能.我已经找到了很多关于如何使用命令的信息,但没有太多关于它实际做了什么(除了一些神秘的评论:http://www.bioconductor.org/help/course-materials/2008/PHSIntro/ lme4Intro-handout-6.pdf).我正在玩以下简单的例子:

library(data.table)
library(lme4)
options(digits=15)

n<-1000
m<-100
data<-data.table(id=sample(1:m,n,replace=T),key="id")
b<-rnorm(m)
data$y<-rand[data$id]+rnorm(n)*0.1
fitted<-lmer(b~(1|id),data=data,verbose=T)
fitted

我理解lmer拟合形式为Y_ {ij} = beta + B_i + epsilon_ {ij}的模型,其中epsilon_ {ij}和B_i分别是具有方差sigma ^ 2和tau ^ 2的独立法线.如果theta = tau/sigma是固定的,我用正确的均值和最小方差计算β的估计值

c = sum_{i,j} alpha_i y_{ij}

哪里

alpha_i = lambda/(1 + theta^2 n_i)
lambda = 1/[\sum_i n_i/(1+theta^2 n_i)]
n_i = number of observations from group i

我还计算了sigma ^ 2的以下无偏估计:

s ^ 2 =\sum_ {i,j} alpha_i(y_ {ij} - c)^ 2 /(1 + theta ^ 2 - lambda)

这些估计似乎与lmer产生的一致.但是,我无法弄清楚在这种情况下如何定义对数似然.我计算了概率密度

pd(Y_{ij}=y_{ij}) = …

我正在尝试拟合混合效果模型,然后使用该模型生成可能具有不同级别的新数据集的估计值.我预计新数据集的估计值将使用估计参数的平均值,但情况似乎并非如此.这是一个最小的工作示例:

library(lme4)
d = data.frame(x = rep(1:10, times = 3),
               y = NA,
               grp = rep(1:3, each = 10))
d$y[d$grp == 1] = 1:10 + rnorm(10)
d$y[d$grp == 2] = 1:10 * 1.5 + rnorm(10)
d$y[d$grp == 3] = 1:10 * 0.5 + rnorm(10)
fit = lmer(y ~ (1+x)|grp, data = d)
newdata = data.frame(x = 1:10, grp = 4)
predict(fit, newdata = newdata, allow.new.levels = TRUE)

在这个例子中,我基本上定义了三组具有不同回归方程(斜率为1,1.5和0.5).但是,当我尝试使用看不见的级别预测新数据集时,我会得到一个恒定的估计值.我本来期望斜率和截距的预期值用于生成这个新数据的预测.我期待错误的事吗？或者,我的代码出了什么问题？

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

当有多个预测变量时,是否可以绘制混合模型的随机截距或斜率？

有了一个预测器,我会这样做:

#generate one response, two predictors and one factor (random effect)
resp<-runif(100,1, 100)
pred1<-c(resp[1:50]+rnorm(50, -10, 10),resp[1:50]+rnorm(50, 20, 5))
pred2<-resp+rnorm(100, -10, 10)
RF1<-gl(2, 50)

#gamm
library(mgcv)
mod<-gamm(resp ~ pred1, random=list(RF1=~1))
plot(pred1, resp, type="n")
for (i in ranef(mod$lme)[[1]]) {
abline(fixef(mod$lme)[1]+i, fixef(mod$lme)[2])
}

#lmer
library(lme4)
mod<-lmer(resp ~ pred1 + (1|RF1))
plot(pred1, resp, type="n")
for (i in ranef(mod)[[1]][,1]) {
abline(fixef(mod)[1]+i, fixef(mod)[2])
}

但是,如果我有这样的模型呢？:

mod<-gamm(resp ~ pred1 + pred2, random=list(RF1=~1))

或者与lmer

mod<-lmer(resp ~ pred1 + pred2 + (1|RF1))

我应该考虑所有系数还是只考虑我正在绘制的变量？

谢谢

我试图将交叉非线性随机效应模型拟合为此问题中提到的线性随机效应模型,并使用该包在此邮件列表中nlme.虽然,无论我尝试什么,我都会收到错误.这是一个例子

library(nlme)

#####
# simulate data
set.seed(18112003)
na <- 30
nb <- 30

sigma_a <- 1
sigma_b <- .5
sigma_res <- .33

n <- na*nb

a <- gl(na,1,n)
b <- gl(nb,na,n)
u <- gl(1,1,n)

x <- runif(n, -3, 3)

y_no_noise <- x + sin(2 * x)
y <- 
  x + sin(2 * x) + 
  rnorm(na, sd = sigma_a)[as.integer(a)] + 
  rnorm(nb, sd = sigma_b)[as.integer(b)] + 
  rnorm(n, sd = sigma_res)

#####
# …

我假设混合效应模型中的随机效应方差对于固定因子的不同水平将有所不同BTyp。

这是我的模特

fm2 <- lme(CA ~ 1 + pF+Tiefe+BTyp+Tiefe:pF+BTyp:pF, data=data2, 
           random = list(~ 1 + pF|Probe))
fm2_Btyphet<-update(fm2, weights=varIdent(form=~1|BTyp))

我Btyp使用lmer函数设法将特定的方差纳入随机效应中，但是此函数不允许考虑组内误差的方差异质性（在我的情况下最好考虑）。我的问题是如何使用lme函数将特定于“ Btyp”的差异纳入随机效应？

在下面您可以看到它如何与lmer功能一起使用。

CA ~ 1 + pF + Tiefe + BTyp + Tiefe:pF + BTyp:pF + 
     (0 + Pind + pF | Probe) + (0 + Bind + pF | Probe) + (0 + Tind + pF | Probe) 


 Data: data2 

 AIC   BIC logLik deviance REMLdev

   21987 22092 -10975    21979 …

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

我对Statsmodels Mixedlm的输出感到有些困惑，希望有人能解释一下。

我有一个庞大的单户住宅数据集，包括每个物业的前两个销售价格/销售日期。我已经对整个数据集进行了地理编码，并获取了每个属性的海拔高度。我试图了解不同城市之间海拔和房地产价格升值之间的关系如何变化。

我使用statsmodels混合线性模型来使价格升值回归高程，同时将其他一些因素保持不变，而城市是我的组类别。

md = smf.mixedlm('price_relative_ind~Elevation+YearBuilt+Sale_Amount_1+LivingSqFt',data=Miami_SF,groups=Miami_SF['City'])

mdf = md.fit()

mdf.random_effects

输入mdf.random_effects将返回系数列表。我能否将此列表解释为实质上每个城市的斜率（即，将高程与销售价格升值相关的单个回归系数）？还是这些结果是每个城市的截距？