1

我正在尝试学习 R 中的分层模型，并且我已经为自己生成了一些示例数据。我在编写多级回归问题的正确语法时遇到问题。

我在商学院生成了一些工资数据。我使薪水与工作年限和教员的出版物总数成线性关系。教师在各个部门，我为每个部门制定了不同的基本工资（截距），并且每个部门的年度加薪（斜率）也不同。这样，我的工资的截距（基本工资）和斜率（多年经验）取决于嵌套级别（部门）和斜率 wrt 另一个解释变量（出版物）不依赖于嵌套级别。在 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)

边界（奇异）拟合：参见 ?isSingular

我想看看部门对工资截距和每年加息的估计，以及作为独立（汇总）发布的影响的估计。现在我根本没有让代码工作。

我知道基本工资和部门的年度加息以及出版物的影响（因为我生成了它）。

Dept    base    Slope
Fin     150000  3000 
Mark    100000  6000 
Ops     120000  4000
IT      140000  5000
Strat   150000  2000

每发表一篇文章，工资就增加一万。

答案： 感谢@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

这些是对原始值的相当好的估计。现在我需要学会评估它们“有多好”。:)

4

1 回答 1

1

(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)

(2)

~ 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 FAQ 有更多关于模型规范的信息

于 2019-08-16T00:46:24.887 回答