Оценка коэффициентов регрессии на уровне смеси c без использования пакета lme4 в R - PullRequest
Новая
       123

Оценка коэффициентов регрессии на уровне смеси c без использования пакета lme4 в R

0 голосов
/ 01 мая 2020

У меня есть двухуровневый набор данных из 37000 экземпляров, который представляет выбор из 199 предметов. Я должен оценить коэффициенты в логистической регрессии c для каждого из 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

1 Ответ

0 голосов
/ 02 мая 2020

Вы можете использовать функцию by() из базового пакета, краткое описание которого "Применить функцию к разделению фрейма данных по факторам" (ref: help(by))

Вот пример использования вашей терминологии для фрейма данных и имен переменных идентификатора субъекта: 

# 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

Если мы go немного дальше, мы также можем выполнить проверку работоспособности , на которой основана каждая модель 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

в обоих случаях.

...