在不使用R中的lme4包的情况下估计混合级逻辑回归系数。

Question

我有一个包含37000个实例的两级数据集，它代表了199个被试的选择。我必须在逻辑回归中估计199个个体的系数。我已经通过子集手动完成了199次，但我想知道是否有更有效的方法可以在不使用lme4包的情况下通过循环获得系数。另外，我应该把系数作为每个科目的变量来计算。

这是我的代码。

### Split of the dataset in each subject ID
mylist <- split(df_merged2, df_merged2$sjind)

### Indication of subject 1 in the first subsetting
df1 <- mylist[[1]]

### Logistic regression

glm1 <- glm(rep ~ reward_v.2 + trans_v.2 + reward_transition, data = df1)   

### Extracting the coefficients
reward_transition <- coef(glm1)[4] 

reward <- coef(glm1)[2] 

transition <- coef(glm1)[3] 

reward<- as.numeric(reward)

reward_transition <- as.numeric(reward_transition)

transition <- as.numeric(transition)

omega <- reward_transition - reward

### Computing the constant coefficients as variables

df1$rewardmix <- 1

df1$rewardmix <- reward

df1$omega <- 1

df1$omega <- omega
df1$transmix <- 1

df1$transmix <- transition

df1$reward_transitionmix <- reward_transition

Answer 1

你可以使用 by() 函数，它的简短描述是 "将函数应用于按因素划分的数据框架" (参考文献) help(by))

下面是一个使用你的术语来说明数据框架和主题ID变量名称的例子。

# Make the simulated data reproducible
set.seed(1717)

# The IDs can be sorted in any order
ids = c('A','B','B','A','A','B','B','B','C','C','C','B','C')

# Sample data frame with: subject ID, target variable (y), input variable (x)
df_merged2 = data.frame(sjind=ids,
                        y=rnorm(length(ids)),
                        x=rnorm(length(ids)))
head(df_merged2)

数据的前6行是这样的：

  sjind          y          x
1     A -1.4548934  1.1004932
2     B -1.7084245 -0.7731208
3     B  2.1004557 -1.6229203
4     A -1.0283021  0.4233806
5     A  0.4133888  1.2398577
6     B -1.4104637  0.3746706

现在使用 by() 函数来拟合每个组的GLM模型。sjind 唯一的值。

glm_by_sjind = by(df_merged2, as.factor(df_merged2$sjind),
                              function(df) glm(y ~ x, data=df))

输出对象 glm_by_sjind 是一个具有以下属性的列表。

• 它的元素数量与列表中的唯一值数量一样多。sjind (本例中为3)
• 它的索引是由以下的唯一值组成的 sjind 变量 "A", "B", "C")
• 每个元素都包含来自 glm() 在输入数据帧的每个分割点上运行（其中分割点由 sjind 唯一值)

因此，例如，你可以要求对主题的回归输出进行总结。"B" 如下所示。

> summary(glm_by_sjind[["B"]])

Call:
glm(formula = y ~ x, data = df)

Deviance Residuals: 
       2         3         6         7         8        12  
-1.40226   1.59040  -0.00186   0.06400  -1.93118   1.68091  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -1.0487     0.7472  -1.404    0.233
x            -0.9605     0.9170  -1.047    0.354

(Dispersion parameter for gaussian family taken to be 2.763681)

    Null deviance: 14.087  on 5  degrees of freedom
Residual deviance: 11.055  on 4  degrees of freedom
AIC: 26.694

Number of Fisher Scoring iterations: 2

---

如果我们再往前走一点，我们还可以做一个 安全检查 每个GLM模型都是基于预期的案例数（即每个模型中的案例数应该等于频率分布的频率）。sjind 变量在输入数据帧中）。)

freq_sjind_in_data = as.list( table(df_merged2$sjind) )
ncases_in_each_glm = lapply( glm_results, function(glm) NROW(glm$data) ) 
all.equal( freq_sjind_in_data,
           ncases_in_each_glm )

其中返回 TRUE.

或者也可以目测一下。

as.data.frame(freq_sjind_in_data)
as.data.frame(ncases_in_each_glm)

其中返回

  A B C
1 3 6 4

在这两种情况下。