Svyglm в обзоре пакетов в R не возвращает стандартные ошибки

Я был бы очень признателен за помощь в этом. Я хотел бы оценить коэффициенты и 95% ДИ для glm, который применяется к обследованию домохозяйств с 2 уровнями (определяемыми дд и чч.число1). Я только недавно наткнулся на пакет опрос.

Я следил за примерами в виньетке для 1) настройки набора данных для учета методов выборки - с помощью svydesign 2) настройки glm с помощью команды svyglm. Для примеров наборов данных:

library(survey)data(api)head(apiclus1)dclus1 <- svydesign(id = ~dnum, weights = ~pw, data = apiclus1)logitmodel <-svyglm(I(sch.wide=="Yes")~awards+comp.imp+enroll+target+hsg+pct.resp+mobility+ell+meals, design=dclus1, family=quasibinomial())summary(logitmodel)

Добавление большого количества переменных выглядит нормально, поэтому я уверен, что пакет работает с хорошим набором данных.

Когда я делаю то же самое со своим набором данных, стандартные ошибки возвращаются с «Inf», если добавляются 3 или 4 переменные, и я не могу понять, почему. Кажется, что это более распространено с факторами. Мне жаль, что я не смог воспроизвести ошибку с другими примерами, но набор данных может быть скачивается здесь.

Итак, используя этот набор данных:

load("balo2_7March17.Rdat")  
dclus1 <- svydesign(id=~dd+hh.num1, weights=~chweight, data = balo2)  
glm1 <- svyglm(out.penta ~ factor(MN18c) + windex5 + age.y, 
          design=dclus1, family=quasibinomial())  
summary(glm1)

Если MN18c является числовым, то выдаются стандартные ошибки, если это множитель (а так и должно быть), стандартными ошибками являются Inf. Не зная, что еще делать, мне нужно попробовать анализ в STATA. Я видел некоторые комментарии о том, что ошибки могут возникнуть при применении к «плохому» набору данных, но что включает в себя «плохой»?


r survey standard-error
person Kath O'Reilly    schedule 09.03.2017    source источник
comment
Привет @Kath O'Reilly, ты смогла решить эту проблему?   -  person student    schedule 05.04.2017

Ответы (1)


arrow_upward
0
arrow_downward

Проблема в том, что в вашей модели нет остаточных степеней свободы. Остаточная df — это плановая df (количество ПЕВ минус количество слоев) за вычетом количества предикторов, которое может легко стать отрицательным, если у вас есть два больших кластера на слой. Это определение остаточного df, вероятно, консервативно, но это не простой вопрос.

>  degf(dclus1)
[1] 5
> glm1$df.resid
[1] 0

Вы можете извлечь стандартные ошибки с помощью

> SE(glm1)
   (Intercept) factor(MN18c)2 factor(MN18c)3 factor(MN18c)4        windex5 
     0.5461374      0.4655331      0.2805168      0.3718879      0.1376936 
         age.y 
     0.1638210

и если вы хотите использовать другие остаточные степени свободы, вы можете указать это summary и получить $p$-значения. В частности, если ни одна из ваших ковариат не находится на уровне кластера, есть разумный аргумент, что регрессия не использует степени свободы, и поэтому для одного параметра за раз вы можете сделать

> summary(glm1, df=degf(dclus1))

Call:
svyglm(formula = out.penta ~ factor(MN18c) + windex5 + age.y, 
    design = dclus1, family = quasibinomial())

Survey design:
svydesign(id = ~dd + hh.num1, weights = ~chweight, data = balo2)

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)     -3.0848     0.5461  -5.648  0.00241 **
factor(MN18c)2  -0.1183     0.4655  -0.254  0.80957   
factor(MN18c)3  -0.4908     0.2805  -1.750  0.14059   
factor(MN18c)4  -0.6137     0.3719  -1.650  0.15981   
windex5          0.2556     0.1377   1.856  0.12256   
age.y            0.9934     0.1638   6.064  0.00176 **

Объединение параметров (например, для проверки трех коэффициентов, составляющих MN18c) более проблематично, и я думаю, что вам как минимум нужно df=degf(clus1)-3+1.

В следующей версии 4.1 пакет будет сообщать о стандартных ошибках в этой ситуации (но не о значениях $p$, если не указано другое df=)

person Thomas Lumley    schedule 06.07.2020