Подходящая отрицательная биноминальная модель

Question

Я пытаюсь запустить отрицательную биномиальную модель для следующей модели.Воспроизведение результатов таблицы 3 в этом документе Связь между реформами закона о вооружениях

di = ln (ni) + β00 + β10 yeari + ei, i = 1997,…, 2013

где, d = общий коэффициент смертности, n = человеко-годы в группе риска и годы с 1997 по 2013 год.

Приведенная выше модель используется для оценки тенденции (измеряется как среднее относительное изменение в смертипоказатели за год) смертности в периоды (1997-2013 гг.)

Как я могу подогнать эту модель в R. Любые предложения по синтаксису в R. Также может кто-нибудь объяснить, что такое β00 и β10?Как мне выбрать эти значения?

Ниже приведены значения в моем фрейме данных

'data.frame':   17 obs. of  13 variables:
 $ X                : int  1 2 3 4 5 6 7 8 9 10 ...
 $ year             : int  1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 .
 $ personyearsatrisk: int  18423037 18607584 18812264 19028802 19274701 19495210 19720737 19932722 20176844 20450966 ...
 $ suicidegun       : int  296 248 261 222 262 213 188 171 145 192 ...
 $ homocidegun      : int  74 56 51 60 48 50 46 18 20 36 ...
 $ massgun          : int  0 0 0 0 0 0 0 0 0 0 ...
 $ suicidenotgun    : int  2351 2393 2229 2170 2208 2115 1976 1948 1951 1952 ...
 $ homocidenotgun   : int  253 243 247 246 268 243 228 147 171 210 ...
 $ suicidetotal     : int  2647 2641 2490 2392 2470 2328 2164 2119 2096 2144 ...
 $ homocidetotal    : int  327 299 298 306 316 293 274 165 191 246 ...
 $ avgsuicidegun    : num  1.61 1.33 1.39 1.17 1.36 ...
 $ avghomocidegun   : num  0.402 0.301 0.271 0.315 0.249 ...
 $ avgsum           : num  2.01 1.63 1.66 1.48 1.61 ...

Ниже приведены мои значения dput (датафрейм)

structure(list(X = 1:17, year = 1997:2013, personyearsatrisk = c(18423037L, 
18607584L, 18812264L, 19028802L, 19274701L, 19495210L, 19720737L, 
19932722L, 20176844L, 20450966L, 20827622L, 21249199L, 21691653L, 
22031750L, 22340024L, 22728254L, 23117353L), suicidegun = c(296L, 
248L, 261L, 222L, 262L, 213L, 188L, 171L, 145L, 192L, 188L, 184L, 
165L, 174L, 145L, 176L, 166L), homocidegun = c(74L, 56L, 51L, 
60L, 48L, 50L, 46L, 18L, 20L, 36L, 29L, 26L, 37L, 39L, 32L, 41L, 
35L), massgun = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L), suicidenotgun = c(2351L, 2393L, 2229L, 
2170L, 2208L, 2115L, 1976L, 1948L, 1951L, 1952L, 2072L, 2141L, 
2196L, 2239L, 2329L, 2374L, 2454L), homocidenotgun = c(253L, 
243L, 247L, 246L, 268L, 243L, 228L, 147L, 171L, 210L, 199L, 229L, 
232L, 202L, 240L, 221L, 184L), suicidetotal = c(2647L, 2641L, 
2490L, 2392L, 2470L, 2328L, 2164L, 2119L, 2096L, 2144L, 2260L, 
2325L, 2361L, 2413L, 2474L, 2550L, 2620L), homocidetotal = c(327L, 
299L, 298L, 306L, 316L, 293L, 274L, 165L, 191L, 246L, 228L, 255L, 
269L, 241L, 272L, 262L, 219L)), .Names = c("X", "year", "personyearsatrisk", 
"suicidegun", "homocidegun", "massgun", "suicidenotgun", "homocidenotgun", 
"suicidetotal", "homocidetotal"), class = "data.frame", row.names = c(NA, 
-17L))

Я рассчитал avgsuicidegun для 100 000 человек в год в основном(suicidegun * 100000) / человеко-годовой риск и аналогично avghomocidegun.

Avgsum - это сумма avgsuicidegun и avghomocidegun

На рисунке ниже показан средний показатель смертности от огнестрельного оружия в год.Пожалуйста, обратитесь к линии тренда с 1997 по 2013 год на этом рисунке:

.

Я пытаюсь найти оценку отношения трендов в годовом коэффициенте смертности на 100000 населения (95% ДИ) за смерть с применением огнестрельного оружия

Answer 1

Я согласен с @Grubbmeister, что ваше описание модели выглядит странно.Вероятно, должно быть:

Для соответствия этой модели:

## compute total gun deaths per year
gundata <- transform(gundata, guntotal=suicidegun+homocidegun)
library(MASS)
g1 <- glm.nb(guntotal ~ 1 + year + offset(log(personyearsatrisk)),
       data=gundata)

Чтобы извлечь нужные коэффициенты:

coef(g1)  
coef(summary(g1))
1-exp(coef(g1)["year"])  ## 0.04903379

Коэффициент "год" равен -0,05027, что означает, что смертность от оружия уменьшается примерно на 5% в год;более точное значение определяется приведенным выше расчетом (т. е. уменьшение составляет 4,9% в год).

Я только кратко ознакомился с бумагой, которую вы связали, но коэффициент здесь (exp(coef(g1)["year"])), по-видимому, согласуется сзначение 0,951, которое они указали для тенденции в период 1997-2013 гг. (это ожидаемое мультипликативное снижение смертности от оружия в течение одного года; отрицательная биномиальная модель использует логарифмическую функцию связи). * * тысяча двадцать-одна

Answer 2

Я не совсем понимаю, как можно использовать отрицательное биномиальное распределение для предиктора, который был преобразован в лог, так как он предназначен для целочисленных данных.

В любом случае ... мы будем работать с тем, что имеем.Я ничего не преобразовывал в журнал, но для этого поместите log() вокруг переменной, которую нужно преобразовать.

Вы должны были также предоставить следующую строку кода:

death$avgsum<-((death$homocidegun+death$suicidegun)/death$personyearsatrisk)*100000

    install.packages("MASS")
    library("MASS")
    mod<-glm.nb(avgsum~year, data= death)
    #to check the residuals
    plot(mod)



#However, I don't think a negative binomial distribution works with this data, so I'd just use a simple linear regression instead:

mod1<-lm(avgsum~year, data= death)

Бета, о которых вы спрашиваете, - это константы в окончательной модели, которые вы можете получить:

summary(mod)

Чтобы получить значения R2 для отрицательной биномиальной модели:

install.packages("jtools")
library(jtools)
summ(mod)