用多层次模型估计相关参数

机器算法验证 r 多层次分析 lme4-nlme
2022-03-31 09:26:12

我想在Stata或R(使用lmer)中估计一个多级模型,其中所有观察的第一级系数都相同,但观察内的系数是相关的。

一个例子看起来像这样:

Yi=β1x1i+β2x2i+β3x3i+...+ε0i
β1=γ1z1+γ2z2+ε1
β2=γ1z1+γ3z3+ε2
β3=γ2z2+γ3z3+ε3
等等,每个 beta 都有方程。

显然,我会为的...做一个分布假设,例如εεN(0,σ2)

x 变量因观测而异,但 z 变量在观测之间没有变化。因此,对于所有观察,参数也是相同的。 γβ

这与我见过的大多数层次模型不同,因为参数在观察中是相关的,而不是取决于观察级别的特征。

作为一个具体的应用,考虑一个模型,其中因变量是学生的考试成绩。x 变量是以前类别的性能度量,变量是这些类别的特征。学生上过相同的课程,但每个班级的学生很少,所以我想汇总系数的估计。因为这些类具有相似的特征,参数可能比参数少得多,并且对那些较低级别的类特征的估计可能会比没有第二级模型的估计 Yzβγββ

同时,我想估计参数,所以代入和估计 y 作为和 x 的函数只会让我走到一半。βγ

估计此类模型的最佳方法是什么?我通常使用 R、Stata 和 Python 进行编程。

3个回答

您是否尝试过使用 Bugs 或 Jags,从 R 中调用其中之一?您似乎正在估计的模型是一个简单的变化斜率模型,预测变量位于第二级。

我会将您的模型重写为:

学生和个班级。假设您的数据是学生类的形式(即重复测量),那么您的模型是:i=1,...nk=1,...K

yiN(β[k]x1,i+δ[k]x2,i+...,σ2)

β[k]N(γ1Z1,k,σβ12) ...

使用 Bugs 或 jags 很容易估计此模型,您可以使用函数 rjags 或 bugs 调用它们。它们在 R2jags 包中,是一个在 R 上拟合多级模型(带有 winBugs)的简单示例。

写出似然函数并最大化怎么样?

这与正常变化系数模型相比有何优势,例如:

fit<-lmer(score~1+vector of class_attributes+vector of student attributes
+(1+vector of class attributes+vector of student attributes)
+(1+vector of student attributes|class)
+(1+vector of class attributes|student))

?

在这个例子中,有一个整体的截距和属性效应,但是每个类都有不同的可能系数,可以通过键入 ranef(fit) 来查看

贝茨关于 lme4 的书的第 3.2 节似乎与您的情况完全相似。

https://r-forge.r-project.org/scm/viewvc.php/*checkout*/www/lMMwR/lrgprt.pdf?revision=656&root=lme4&pathrev=656

更新(我更新了上面的代码行):

我还运行了这些行来尝试模拟您的情况,但没有任何学生属性

library(lme4)
n<-100 #class size
pool<-200 #student pool size
class=c(rep(1,n), rep(2,n), rep(3,n))
min_in_class=c(rep(45,n), rep(60,n), rep(90,n))
min_hw=c(rep(90,n), rep(60,n), rep(60,n))
student_id=c(sample(1:pool,n), sample(1:pool,n), sample(1:pool,n))
performance=55+10*class +.1*min_in_class    +.2*min_hw+ -.001*min_in_class*min_hw   +rnorm(3*n, 0,10)
df<-data.frame(class=as.factor(class), min_in_class, min_hw, student_id=as.factor(student_id), performance)
library(reshape2)
melted<-melt(df, id.vars=c('student_id', 'class'))
casted<-dcast(melted, student_id~class+variable)
casted$score<-rowMeans(casted[,c(4,7,10)],na.rm=T)+rnorm(nrow(casted),0,5)
df$score<-casted$score[match(df$student_id, casted$student_id)]

我认为您需要尝试做的是:

fit<-lmer(score~1+min_in_class+min_hw+(1|class)+(1+min_in_class+min_hw|student_id), data=df)

我用各种班级规模和池子运行它,但没有得到我期望的结果;但也许有更多的课程,事情会看起来更好。