我正在尝试学习R中的分层模型,并为自己生成了一些示例数据。我在编码多级回归问题的正确语法方面遇到问题。
我生成了一些商学院的工资数据。我使工资线性取决于教师的就业年限和出版物总数。教师在各个部门,我使每个部门的基本工资(截距(不同,每个部门的年度加息(斜坡(也不同。这样,我就有了工资的截距(基本工资(和斜率(w.r.t年数经验(取决于嵌套级别(部门(和斜率w.r.t另一个解释变量(出版物(不依赖于嵌套级别。在 R 中对此进行建模的正确语法是什么?
这是我的数据
Data <-data.frame(Sl_No = c(1:40),
+ Dept = as.factor(sample(c("Mark","IT","Fin"),40,replace = TRUE)),
+ Years = round(runif(40,1,10)))
pubs <-round(Data$Years*runif(40,1,3))
Data$Pubs <- pubs
lookup_table<-data.frame(Dept = c("Mark","IT","Fin","Strat","Ops"),
+ base = c(100000,140000,150000,150000,120000),
+ slope = c(6000,5000,3000,2000,4000))
Data <- merge(Data,lookup_table,by = 'Dept')
salary <-Data$base+Data$slope*Data$Years+Data$Pubs*10000+rnorm(length(Data$Dept))*10000
Data$base<-NULL
Data$slope<-NULL
我尝试了以下方法:
1(
multilevel_model<-lmer(Salary~1|Dept+Pubs+Years|Dept, data = Data)
model.matrix.default(eval(substitute(~foo, list(foo = x[[2]](((, : 模型框架和公式在 model.matrix(( 中不匹配
2(
multilevel_model<-lmer(`Salary`~ Dept + `Pubs`+`Years`|Dept , data = Data)
边界(奇异(拟合:参见 ?是单数
我想看看部门对工资截获和年度加薪的估计,以及对独立(合并(出版效果的估计。现在我根本没有让代码工作。
我知道基本工资和部门每年的涨幅以及出版物的效果(因为我产生了它(。
Dept base Slope
Fin 150000 3000
Mark 100000 6000
Ops 120000 4000
IT 140000 5000
Strat 150000 2000
每出版一份出版物,工资就会增加10,000。
答:感谢@Ben在这里的回答,我认为正确的模型是
multilevel_model<-lmer(Salary~(1|Dept)+ Pubs +(0+Years|Dept), data = Data)
这通过运行为我提供了以下固定效果
summary(multilevel_model)
Fixed effects:
Estimate Std. Error t value
(Intercept) 131667.4 10461.0 12.59
Pubs 10235.0 550.8 18.58
Correlation of Fixed Effects:
Pubs -0.081
部门级系数如下:
coef(multilevel_model)
$Dept
Years (Intercept) Pubs
Fin 3072.5133 148757.6 10235.02
IT 5156.6774 136710.7 10235.02
Mark 5435.8301 102858.3 10235.02
Ops 3433.1433 118287.1 10235.02
Strat 963.9366 151723.1 10235.02
这些是对原始值的相当不错的估计。现在我需要学会评估他们"有多好"。:)
(1(
multilevel_model<-lmer(`Total Salary`~ 1|Dept +
`Publications`+`Years of Exp`|Dept , data = sample_data)
我无法立即诊断为什么这会给出语法错误,但通常建议在随机效应项周围使用括号,因为|
运算符在公式中具有很高的优先级。 因此,响应/右侧 (RHS( 公式
~ (1|Dept) + (`Publications`+`Years of Exp`|Dept)
可能会起作用,但会有问题,因为两个术语都包含相同的截距项:如果你想这样做,你可能需要
~ (1|Dept) + (0+`Publications`+`Years of Exp`|Dept)
(二(
~ Dept + `Publications`+`Years of Exp`|Dept
将相同的变量(Dept
(放在柱线的左侧和右侧没有任何意义。
您可能应该使用
~ pubs + years_exp + (1 + years_exp|Dept)
由于原则上出版的效果可能因部门而异,因此最大模型将是
~ pubs + years_exp + (1 + pubs + years_exp|Dept)
- 包含没有相应固定效应的随机效应很少有意义。
- 请注意,即使您拥有正确的模型,也可能获得单数拟合;请参见
?isSingular
手册页。 - 如果上面列出的 18 个观测值代表您的整个数据集,则很可能太小而无法成功拟合最大模型。经验法则是,每个估计参数需要 10-20 个观测值,最大模型具有 (截距 + 2 个固定效应参数 + (3*4(/2=6 个随机效应参数(= 9 个参数。(由于它是模拟的,因此您可以轻松模拟大数据集...
- 我建议重命名数据框中的变量,这样您就不必为带有空格的反引号保护变量名称大惊小怪......
GLMM常见问题解答有更多关于型号规格的信息