تجزیه و تحلیل رگرسیون پواسن در R

در مدل متغیرهای تعداد دفعات مشاهده شده رگرسیون پواسون استفاده می شود . برای وارد کردن کد های این قسمت باید پکیج هایی بر روی نرم افزار R نصب شود که این پکیج ها عبارتند از :

require(ggplot2) require(sandwich) require(msm)

  نمونه هایی از رگرسیون پواسون مثال 1. تعداد افرادی که توسط قاطر یا اسب لگد در ارتش پروس در سال کشته شدند. . Ladislaus. Bortkiewicz اطلاعات از 20 جلد از Preussischen Statistik جمع آوری شده. این داده ها در 10 سپاه از ارتش پروس در اواخر 1800s در طول 20 سال جمع آوری شد. مثال 2. تعداد افراد در خط در مقابل شما در فروشگاه مواد غذایی. پیش بینی ممکن است تعدادی از اقلام در حال حاضر در قیمت با تخفیف ویژه و اینکه آیا یک رویداد خاص (به عنوان مثال، یک روز تعطیل، یک رویداد ورزشی بزرگ) است سه یا کمتر روز دور ارائه شده است. مثال 3 تعداد جوایز کسب شده توسط دانش آموزان در یک دبیرستان. عوامل تعیین کننده تعداد جوایز به دست آورده شامل نوع برنامه است که در آن دانش آموزان ثبت نام کرده بود (به عنوان مثال، حرفه ای، عمومی و یا آکادمیک) و نمره در امتحان نهایی خود را در ریاضی. شرح داده به منظور تصویر، ما شبیه سازی یک مجموعه داده برای مثال 3 کنید. در این مثال، num_awards متغیر نتیجه است و نشان می دهد که تعدادی از جوایز کسب شده توسط دانش آموزان در یک دبیرستان در یک سال، math یک متغیر پیش بینی کننده مستمر است و نمرات دانش آموزان در امتحان ریاضی نهایی خود را، و prog یک متغیر پیش بینی کننده قطعی است با سه سطح نشان می دهد نوع برنامه است که در آن دانش آموزان انجام شد. آن را به عنوان 1 = "عمومی"، 2 = "علمی" و 3 = "حرفه ای" کد شده. اجازه دهید با بارگذاری داده ها و به دنبال برخی از آمار توصیفی شروع می شود.

p <- read.csv("http://www.ats.ucla.edu/stat/data/poisson_sim.csv") p <- within(p, { prog <- factor(prog, levels = 1:3, labels = c("General", "Academic", "Vocational")) id <- factor(id) }) summary(p)

## id num_awards prog math ## 1 : 1 Min. :0.00 General : 45 Min. :33.0 ## 2 : 1 1st Qu.:0.00 Academic :105 1st Qu.:45.0 ## 3 : 1 Median :0.00 Vocational: 50 Median :52.0 ## 4 : 1 Mean :0.63 Mean :52.6 ## 5 : 1 3rd Qu.:1.00 3rd Qu.:59.0 ## 6 : 1 Max. :6.00 Max. :75.0 ## (Other):194

  هر متغیر دارای 200 مشاهدات معتبر و توزیع آنها به نظر می رسد کاملا معقول. میانگین بی قید و شرط و واریانس متغیر نتیجه ما بسیار متفاوت است. مدل ما فرض می شود که این ارزش ها، مشروط بر متغیرهای پیش بینی خواهد شد، برابر (یا حداقل تقریبا تا) باشد. ما می توانیم با استفاده از tapply تابع برای نمایش خلاصه آمارها بر اساس نوع برنامه است. جدول زیر متوسط تعداد جوایز اساس نوع برنامه را نشان می دهد و به نظر می رسد نشان می دهد که نوع برنامه کاندید خوبی برای پیش بینی تعدادی از جوایز، متغیر نتیجه ما است، چرا که مقدار متوسط از نتیجه به نظر می رسد توسط متفاوت prog . علاوه بر این، ابزار و واریانس در هر سطح prog - به معنای مشروط و واریانس - مشابه هستند. هیستوگرام مشروط توسط نوع برنامه از هم جدا شده رسم شده است برای نشان دادن توزیع است.

with(p, tapply(num_awards, prog, function(x) { sprintf("M (SD) = %1.2f (%1.2f)", mean(x), sd(x)) }))

## General Academic Vocational ## "M (SD) = 0.20 (0.40)" "M (SD) = 1.00 (1.28)" "M (SD) = 0.24 (0.52)" ggplot(p, aes(num_awards, fill = prog)) + geom_histogram(binwidth = 0.5, position = "dodge")

  روشهای تجزیه و تحلیل در زیر لیستی از برخی از روش های تجزیه و تحلیل ممکن است شما مواجه می شوند می باشد. برخی از روش های ذکر شده کاملا معقول است، در حالی که دیگران کرده اند یا از نفع کاهش یافته است و یا دارای محدودیت. • رگرسیون پواسون - رگرسیون پواسون است که اغلب برای اطلاعات تعداد مدل سازی استفاده می شود. رگرسیون پواسون تا به تعدادی از الحاقات برای مدل تعداد دفعات مشاهده شده مفید است. • رگرسیون دو جمله ای منفی - رگرسیون دو جمله ای منفی می تواند برای بیش از پراکنده داده های شمارش مورد استفاده قرار گیرد، این است که زمانی که واریانس شرطی بیش از میانگین شرطی. این را می توان به عنوان یک تعمیم رگرسیون پواسون در نظر گرفته از آن است که ساختار متوسط همان رگرسیون پواسون و آن را تا یک پارامتر اضافی به مدل بیش از پراکندگی. اگر توزیع مشروط متغیر نتیجه بیش از حد پراکنده، فاصله اطمینان برای رگرسیون دو جمله ای منفی به نظر می رسد باریک تر در مقایسه با کسانی که از رگرسیون پواسون است . • مدل رگرسیون صفر متورم - مدل صفر متورم به حساب تلاش برای صفر اضافی می باشد. به عبارت دیگر، دو نوع از صفر تصور می شود در داده ها، "صفر واقعی" و "صفر اضافی" وجود داشته باشد. مدل های صفر متورم برآورد دو معادله به طور همزمان، یکی برای مدل تعداد دفعات مشاهده شده و یکی برای صفر اضافی می باشد. • رگرسیون OLS - تعداد متغیرهای گاهی اوقات ورود تبدیل و با استفاده از رگرسیون OLS استفاده شد. بسیاری از مسائل را با این روش بوجود می آیند، از جمله از دست دادن داده ها با توجه به ارزش های تعریف نشده تولید با در نظر گرفتن ورود به سیستم از صفر (که تعریف نشده) و تخمین های مغرضانه است . • رگرسیون پواسون در این مرحله، ما آماده هستیم برای انجام تجزیه و تحلیل مدل پواسون با ما می glm عملکرد. ما متناسب با مدل و ذخیره آن در جسم m1 و یک خلاصه ای از مدل در همان زمان است .

summary(m1 <- glm(num_awards ~ prog + math, family = "poisson", data = p)) ## ## Call: ## glm(formula = num_awards ~ prog + math, family = "poisson", data = p) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.204 -0.844 -0.511 0.256 2.680 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -5.2471 0.6585 -7.97 1.6e-15 *** ## progAcademic 1.0839 0.3583 3.03 0.0025 ** ## progVocational 0.3698 0.4411 0.84 0.4018 ## math 0.0702 0.0106 6.62 3.6e-11 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 287.67 on 199 degrees of freedom ## Residual deviance: 189.45 on 196 degrees of freedom ## AIC: 373.5 ## ## Number of Fisher Scoring iterations: 6

  کامرون و Trivedi (2009) با استفاده از خطاهای قوی استاندارد برای برآورد پارامتر برای کنترل نقض خفیف فرض توزیع که واریانس برابر است با متوسط توصیه می شود. ما با استفاده از R بسته بندی sandwich زیر برای به دست آوردن خطا قوی استاندارد و محاسبه مقدار p درآمده است. همراه با ارزش P، ما نیز فاصله اطمینان 95٪ با استفاده از برآورد پارامترها و خطاهای استاندارد قوی خود را محاسبه می شود

cov.m1 <- vcovHC(m1, type = "HC0") std.err <- sqrt(diag(cov.m1)) r.est <- cbind(Estimate = coef(m1), `Robust SE` = std.err, `Pr(>|z|)` = 2 * pnorm(abs(coef(m1)/std.err), lower.tail = FALSE), LL = coef(m1) - 1.96 * std.err, UL = coef(m1) + 1.96 * std.err)

r.est

## Estimate Robust SE Pr(>|z|) LL UL ## (Intercept) -5.24712 0.64600 4.567e-16 -6.5133 -3.98097 ## progAcademic 1.08386 0.32105 7.355e-04 0.4546 1.71311 ## progVocational 0.36981 0.40042 3.557e-01 -0.4150 1.15463 ## math 0.07015 0.01044 1.784e-11 0.0497 0.09061

  • خروجی با انعکاس فراخوانی تابع شروع می شود. اطلاعات در مورد باقیمانده انحراف است در کنار نمایش داده می شود. باقیمانده انحراف تقریبا به طور معمول اگر مدل correctly.In مثال ما مشخص شده توزیع شده، نشان می دهد یک کمی از skeweness از متوسط است نه کاملا صفر است. • بعدی ضرایب رگرسیون پواسون برای هر یک از متغیرهای همراه با انحراف استاندارد، Z-نمرات، مقدار p و فاصله اطمینان 95٪ برای ضرایب آمده است. ضریب برای math 0.07 است. این بدان معنی است که تعداد ورود به سیستم انتظار برای افزایش یک واحد در math 0.07 است. متغیر شاخص progAcademic بین پروژسترون = "علمی" و پروژسترون = "عمومی"، تعداد ورود به سیستم انتظار برای پروژسترون = "علمی" را افزایش می دهد در حدود 1.1 مقایسه می کند. متغیر شاخص prog.Vocational تفاوت انتظار می رود در تعداد دفعات ورود به سیستم (\ (\ حدود 0.37 \)) بین پروژسترون = "حرفه ای" و گروه مرجع(پروژسترون = "عمومی" )است. • اطلاعات در مورد انحراف نیز ارائه شده است. ما می توانیم انحراف باقی مانده به انجام خوبی از آزمون مناسب برای مدل کلی استفاده کنید. انحراف باقی مانده تفاوت بین انحراف از مدل فعلی و حداکثر انحراف از مدل ایده آل که در آن مقادیر پیش بینی شده یکسان به مشاهده می باشد. بنابراین، اگر تفاوت باقی مانده به اندازه کافی کوچک است، به خوبی از آزمون مناسب نمی تواند قابل توجه باشد، نشان می دهد که مدل متناسب با داده ها. ما نتیجه گیری کرد که مدل متناسب به خوبی به این دلیل که آزمون مجذور کای خوبی، از، مناسب است از نظر آماری معنی دار نبود. اگر آزمون آماری معنی دار بوده است، این امر نشان می دهد که داده ها به مدل متناسب نیست و. در این وضعیت، ما ممکن است سعی کنید برای تعیین اینکه آیا وجود دارد حذف متغیرهای پیش بینی، اگر فرض خطی بودن ما را نگه می دارد و / یا اگر موضوع وجود دارد بیش از پراکندگی است.

with(m1, cbind(res.deviance = deviance, df = df.residual, p = pchisq(deviance, df.residual, lower.tail = FALSE)))

## res.deviance df p ## [1,] 189.4 196 0.6182

  ما همچنین می تواند اثر کلی آزمایش prog با مقایسه انحراف از مدل کامل با انحراف از مدل حذف prog . درجه آزادی آزمون دو کای دو نشان می دهد که prog ، گرفته شده با هم، پیش بینی آماری معنی دار است num_awards

## update m1 model dropping prog m2 <- update(m1, . ~ . - prog) ## test model differences with chi square test anova(m2, m1, test = "Chisq") ## Analysis of Deviance Table ## ## Model 1: num_awards ~ math ## Model 2: num_awards ~ prog + math ## Resid. Df Resid. Dev Df Deviance Pr(>Chi) ## 1 198 204 ## 2 196 189 2 14.6 0.00069 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  گاهی اوقات ممکن است بخواهید به ارائه نتایج رگرسیون به نسبت میزان حادثه و خطاهای استاندارد خود، همراه با اعتماد به نفس فاصله است. برای محاسبه خطای استاندارد برای نسبت میزان حادثه، ما به روش دلتا استفاده کنید. برای این منظور، ما را با استفاده از تابع deltamethod اجرا در R بسته msm است .

s <- deltamethod(list(~exp(x1), ~exp(x2), ~exp(x3), ~exp(x4)), coef(m1), cov.m1)

## exponentiate old estimates dropping the p values rexp.est <- exp(r.est[, -3]) ## replace SEs with estimates for exponentiated coefficients rexp.est[, "Robust SE"] <- s

rexp.est ## Estimate Robust SE LL UL ## (Intercept) 0.005263 0.00340 0.001484 0.01867 ## progAcademic 2.956065 0.94904 1.575551 5.54620 ## progVocational 1.447458 0.57959 0.660335 3.17284 ## math 1.072672 0.01119 1.050955 1.09484

خروجی بالا نشان می دهد که این نرخ حادثه برای پروژسترون = "علمی" 2.96 برابر نرخ حادثه برای گروه مرجع (پروژسترون = "عمومی") است. به همین ترتیب، نرخ حادثه برای پروژسترون = "حرفه ای" 1.45 برابر نرخ حادثه برای گروه مرجع برگزاری متغیرهای دیگر در ثابت است. تغییر درصد در میزان حادثه از num_awards برای هر واحد افزایش است 7٪ math . برای کسب اطلاعات بیشتر در معیارهای مختلف است که در آن نتایج را می توان ارائه و تفسیر از جمله، لطفا مدل رگرسیون برای متغیرهای وابسته طبقه با استفاده از STATA، ویرایش دوم توسط J. اسکات لانگ و جرمی Freese (2006) را ببینید. گاهی اوقات، ما ممکن است بخواهید به در وسیله حاشیه ای انتظار می رود نگاه کنید. به عنوان مثال، چه تعداد مورد انتظار برای هر نوع برنامه برگزاری نمره ریاضی در میانگین کلی آن هستند؟ برای پاسخ به این سوال، ما می توانید با استفاده از را predict عملکرد. اول از همه، ما یک داده کوچک مجموعه به درخواست را predict عملکرد به آن است .

(s1 <- data.frame(math = mean(p$math), prog = factor(1:3, levels = 1:3, labels = levels(p$prog)))) ## math prog ## 1 52.65 General ## 2 52.65 Academic ## 3 52.65 Vocational predict(m1, s1, type = "response", se.fit = TRUE) ## $fit ## 1 2 3 ## 0.2114 0.6249 0.3060 ## ## $se.fit ## 1 2 3 ## 0.07050 0.08628 0.08834 ## ## $residual.scale ## [1] 1

در خروجی بالا، ما می بینیم که تعداد پیش بینی رویدادها برای سطح 1 از prog است در مورد 0.21، برگزاری math در متوسط آن است. تعداد پیش بینی رویدادها برای سطح 2 از prog در 0.62 بالاتر است، و تعداد پیش بینی رویدادها برای سطح 3 از prog در حدود 0.31. نسبت از این تعداد پیش بینی شده (\ (\ FRAC {.625} {0.211} = 2.96 \) \ (\ FRAC {.306} {0.211} = 1.45 \)) با آنچه ما دیدیم به دنبال در IRR. است. ما همچنین می توانیم تعداد پیش بینی شده از وقایع با دستورات زیر نمودار. نمودار نشان می دهد که بسیاری از جوایز برای کسانی که در برنامه تحصیلی (پروژسترون = 2) پیش بینی شده است، به خصوص اگر دانشجو نمره ریاضی بالا. پایین ترین تعداد جوایز پیش بینی شده است برای آن دسته از دانش آموزان در برنامه عمومی (پروژسترون=1) نمودار پوشش خطوط از ارزش مورد انتظار را بر روی نقاط واقعی، اگر چه مقدار کمی از سر و صدا به صورت تصادفی به صورت عمودی به کاهش overplotting اضافه شد.

## calculate and store predicted values p$phat <- predict(m1, type = "response")

## order by program and then by math p <- p[with(p, order(prog, math)), ]

## create the plot ggplot(p, aes(x = math, y = phat, colour = prog)) + geom_point(aes(y = num_awards), alpha = 0.5, position = position_jitter(h = 0.2)) + geom_line(size = 1) + labs(x = "Math Score", y = "Expected number of awards")

• هنگامی که به نظر می رسد مسئله پراکندگی وجود دارد، ما ابتدا باید بررسی کنید که آیا مدل ما است به درستی مشخص شده، مانند متغیرهای حذف شده و فرم های عملکردی است. به عنوان مثال، اگر ما پیش بینی متغیر حذف prog در مثال بالا، مدل ما به نظر می رسد مشکل با بیش از پراکندگی داشته باشند. به عبارت دیگر، یک مدل misspecified می تواند یک علامت شبیه به یک مشکل بیش از پراکندگی در حال حاضر. نکات • هنگامی که به نظر می رسد مسئله پراکندگی وجود دارد، ما ابتدا باید بررسی کنید که آیا مدل ما است به درستی مشخص شده، مانند متغیرهای حذف شده و فرم های عملکردی است. به عنوان مثال، اگر ما پیش بینی متغیر حذف prog در مثال بالا، مدل ما به نظر می رسد مشکل با بیش از پراکندگی داشته باشند. به عبارت دیگر، یک مدل misspecified می تواند یک علامت شبیه به یک مشکل بیش از پراکندگی در حال حاضر. • با فرض این که مدل به درستی مشخص شده، این فرض که واریانس شرطی به میانگین شرطی برابر است باید بررسی می شود. تستهای مختلف از جمله آزمون نسبت احتمال بیش از پراکندگی آلفا پارامتر های در حال اجرا به همان مدل با استفاده از توزیع دو جمله ای منفی وجود دارد. R بسته pscl (علوم سیاسی آزمایشگاه محاسباتی، دانشگاه استنفورد) فراهم می کند بسیاری از توابع برای دو جمله ای و تعداد داده ها از جمله odTest برای آزمایش بر روی پراکندگی. • یکی از دلایل شایع بیش از پراکندگی صفر اضافی، که به نوبه خود توسط یک فرایند تولید داده های اضافی تولید شده است. در این وضعیت، مدل صفر متورم باید در نظر گرفته شود. • اگر فرایند تولید داده ها برای هر 0s و (مانند تعداد روز را صرف در بیمارستان) اجازه نمی دهد، پس از آن یک مدل صفر کوتاه ممکن است مناسب تر باشد. • تعداد داده ها اغلب یک متغیر قرار گرفتن در معرض، که نشان می دهد چند بار این رویداد می تواند رخ داده است. این متغیر باید به یک مدل پواسون با استفاده از گنجانیده offset گزینه. • متغیر پیامد در یک رگرسیون پواسون می تواند اعداد منفی را نداشته باشند، و قرار گرفتن در معرض می تواند 0s و ندارد. • رگرسیون پواسون است که از طریق برآورد احتمال حداکثر برآورد شده است. این حالت معمولا نیاز به حجم نمونه بزرگ است.