原文链接：http://tecdat.cn/?p=6560

原文出处：拓端数据部落公众号

读取数据


summary(eba1977)
##          city      age         pop             cases
##  Fredericia:6   40-54:4   Min.   : 509.0   Min.   : 2.000
##  Horsens   :6   55-59:4   1st Qu.: 628.0   1st Qu.: 7.000
##  Kolding   :6   60-64:4   Median : 791.0   Median :10.000
##  Vejle     :6   65-69:4   Mean   :1100.3   Mean   : 9.333
##                 70-74:4   3rd Qu.: 954.8   3rd Qu.:11.000
##                 75+  :4   Max.   :3142.0   Max.   :15.000

普通 Poisson model

glm1 <- glm(formula = cases ~ age + city + offset(log(pop)),family  = poisson(link = "log"),data    = eba1977)
summary(glm1)
##
## Call:
## glm(formula = cases ~ age + city + offset(log(pop)), family = poisson(link = "log"),
##     data = eba1977)
##
## Deviance Residuals:
##      Min        1Q    Median        3Q       Max
## -2.63573  -0.67296  -0.03436   0.37258   1.85267
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -5.6321     0.2003 -28.125  < 2e-16 ***
## age55-59      1.1010     0.2483   4.434 9.23e-06 ***
## age60-64      1.5186     0.2316   6.556 5.53e-11 ***
## age65-69      1.7677     0.2294   7.704 1.31e-14 ***
## age70-74      1.8569     0.2353   7.891 3.00e-15 ***
## age75+        1.4197     0.2503   5.672 1.41e-08 ***
## cityHorsens  -0.3301     0.1815  -1.818   0.0690 .
## cityKolding  -0.3715     0.1878  -1.978   0.0479 *
## cityVejle    -0.2723     0.1879  -1.450   0.1472
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 129.908  on 23  degrees of freedom
## Residual deviance:  23.447  on 15  degrees of freedom
## AIC: 137.84
##
## Number of Fisher Scoring iterations: 5

Stan

数据

## Model matrix
modMat <- as.data.frame(model.matrix(glm1))
modMat$offset <- log(eba1977$pop)
names(modMat) <- c("intercept", "age55_59", "age60_64", "age65_69", "age70_74", "age75plus", "cityHorsens", "cityKolding", "cityVejle", "offset")dat   <- as.list(modMat)
dat$y <- eba1977$cases
dat$N <- nrow(modMat)
dat$p <- ncol(modMat) - 1
## Load Stan file
fileName <- "./poisson.stan"
stan_code <- readChar(fileName, file.info(fileName)$size)
cat(stan_code)
## Run Stan
resStan <- stan(model_code = stan_code, data = dat,chains = 3, iter = 3000, warmup = 500, thin = 10)
##
## TRANSLATING MODEL 'stan_code' FROM Stan CODE TO C++ CODE NOW.
## COMPILING THE C++ CODE FOR MODEL 'stan_code' NOW.
## In file included from file60814bc1cb78.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/model/model_header.hpp:17:
## In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/agrad/rev.hpp:5:
## /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/agrad/rev/chainable.hpp:87:17: warning: 'static' function 'set_zero_all_adjoints' declared in header file should be declared 'static inline' [-Wunneeded-internal-declaration]
##     static void set_zero_all_adjoints() {
##                 ^
## In file included from file60814bc1cb78.cpp:8:
## In file included from /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/model/model_header.hpp:21:
## /Library/Frameworks/R.framework/Versions/3.1/Resources/library/rstan/include//stansrc/stan/io/dump.hpp:26:14: warning: function 'product' is not needed and will not be emitted [-Wunneeded-internal-declaration]
##       size_t product(std::vector<size_t> dims) {
##              ^
## 2 warnings generated.
##
## SAMPLING FOR MODEL 'stan_code' NOW (CHAIN 1).
##
## Iteration:    1 / 3000 [  0%]  (Warmup)
## Iteration:  300 / 3000 [ 10%]  (Warmup)
## Iteration:  501 / 3000 [ 16%]  (Sampling)
## Iteration:  800 / 3000 [ 26%]  (Sampling)
## Iteration: 1100 / 3000 [ 36%]  (Sampling)
## Iteration: 1400 / 3000 [ 46%]  (Sampling)
## Iteration: 1700 / 3000 [ 56%]  (Sampling)
## Iteration: 2000 / 3000 [ 66%]  (Sampling)
## Iteration: 2300 / 3000 [ 76%]  (Sampling)
## Iteration: 2600 / 3000 [ 86%]  (Sampling)
## Iteration: 2900 / 3000 [ 96%]  (Sampling)
## Iteration: 3000 / 3000 [100%]  (Sampling)
## #  Elapsed Time: 0.142295 seconds (Warm-up)
## #                0.543612 seconds (Sampling)
## #                0.685907 seconds (Total)
##
##
## SAMPLING FOR MODEL 'stan_code' NOW (CHAIN 2).
##
## Iteration:    1 / 3000 [  0%]  (Warmup)
## Iteration:  300 / 3000 [ 10%]  (Warmup)
## Iteration:  501 / 3000 [ 16%]  (Sampling)
## Iteration:  800 / 3000 [ 26%]  (Sampling)
## Iteration: 1100 / 3000 [ 36%]  (Sampling)
## Iteration: 1400 / 3000 [ 46%]  (Sampling)
## Iteration: 1700 / 3000 [ 56%]  (Sampling)
## Iteration: 2000 / 3000 [ 66%]  (Sampling)
## Iteration: 2300 / 3000 [ 76%]  (Sampling)
## Iteration: 2600 / 3000 [ 86%]  (Sampling)
## Iteration: 2900 / 3000 [ 96%]  (Sampling)
## Iteration: 3000 / 3000 [100%]  (Sampling)
## #  Elapsed Time: 0.13526 seconds (Warm-up)
## #                0.517139 seconds (Sampling)
## #                0.652399 seconds (Total)
##
##
## SAMPLING FOR MODEL 'stan_code' NOW (CHAIN 3).
##
## Iteration:    1 / 3000 [  0%]  (Warmup)
## Iteration:  300 / 3000 [ 10%]  (Warmup)
## Iteration:  501 / 3000 [ 16%]  (Sampling)
## Iteration:  800 / 3000 [ 26%]  (Sampling)
## Iteration: 1100 / 3000 [ 36%]  (Sampling)
## Iteration: 1400 / 3000 [ 46%]  (Sampling)
## Iteration: 1700 / 3000 [ 56%]  (Sampling)
## Iteration: 2000 / 3000 [ 66%]  (Sampling)
## Iteration: 2300 / 3000 [ 76%]  (Sampling)
## Iteration: 2600 / 3000 [ 86%]  (Sampling)
## Iteration: 2900 / 3000 [ 96%]  (Sampling)
## Iteration: 3000 / 3000 [100%]  (Sampling)
## #  Elapsed Time: 0.120931 seconds (Warm-up)
## #                0.509901 seconds (Sampling)
## #                0.630832 seconds (Total)
## Show traceplot
traceplot(resStan, pars = c("beta"), inc_warmup = TRUE)

比较

## Frequentist
tableone::ShowRegTable(glm1, exp = FALSE)
##             beta [confint]       p
## (Intercept) -5.63 [-6.04, -5.26] <0.001
## age55-59     1.10 [0.61, 1.59]   <0.001
## age60-64     1.52 [1.07, 1.98]   <0.001
## age65-69     1.77 [1.32, 2.22]   <0.001
## age70-74     1.86 [1.40, 2.32]   <0.001
## age75+       1.42 [0.93, 1.91]   <0.001
## cityHorsens -0.33 [-0.69, 0.03]   0.069
## cityKolding -0.37 [-0.74, -0.00]  0.048
## cityVejle   -0.27 [-0.64, 0.09]   0.147
## Bayesian
print(resStan, pars = c("beta"))
## Inference for Stan model: stan_code.
## 3 chains, each with iter=3000; warmup=500; thin=10;
## post-warmup draws per chain=250, total post-warmup draws=750.
##
##          mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
## beta[1] -5.66    0.01 0.21 -6.13 -5.80 -5.64 -5.51 -5.29   655    1
## beta[2]  1.11    0.01 0.25  0.60  0.95  1.11  1.28  1.60   750    1
## beta[3]  1.53    0.01 0.23  1.10  1.38  1.51  1.68  2.00   750    1
## beta[4]  1.77    0.01 0.25  1.30  1.60  1.76  1.94  2.24   750    1
## beta[5]  1.87    0.01 0.24  1.40  1.71  1.86  2.02  2.37   750    1
## beta[6]  1.42    0.01 0.25  0.94  1.25  1.42  1.58  1.95   631    1
## beta[7] -0.33    0.01 0.18 -0.69 -0.45 -0.32 -0.21  0.03   703    1
## beta[8] -0.37    0.01 0.19 -0.74 -0.50 -0.38 -0.24 -0.01   664    1
## beta[9] -0.28    0.01 0.19 -0.66 -0.40 -0.27 -0.15  0.09   698    1
##
## Samples were drawn using NUTS(diag_e) at Mon Apr 13 21:43:02 2015.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).

非常感谢您阅读本文，有任何问题请在下方留言！

拓端tecdat|R语言stan泊松回归Poisson regression相关推荐

  1. 拓端tecdat|R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险

    最近我们被客户要求撰写关于冠心病风险的研究报告,包括一些图形和统计输出. 相关视频:R语言逻辑回归(Logistic回归)模型分类预测病人冠心病风险 逻辑回归Logistic模型原理和R语言分类预测冠 ...

  2. 拓端tecdat|R语言用LOESS(局部加权回归)季节趋势分解（STL）进行时间序列异常检测

    最近我们被客户要求撰写关于LOESS(局部加权回归)的研究报告,包括一些图形和统计输出. 这篇文章描述了一种对涉及季节性和趋势成分的时间序列的中点进行建模的方法.我们将对一种叫做STL的算法进行研究, ...

  3. 拓端tecdat|R语言线性回归和时间序列分析北京房价影响因素可视化案例

    最近我们被客户要求撰写关于北京房价影响因素的研究报告,包括一些图形和统计输出. 目的 房价有关的数据可能反映了中国近年来的变化: 人们得到更多的资源(薪水),期望有更好的房子 人口众多 独生子女政策: ...

  4. 拓端tecdat|R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系

    最近我们被客户要求撰写关于向量误差修正模型的研究报告,包括一些图形和统计输出. 向量自回归模型估计的先决条件之一是被分析的时间序列是平稳的.但是,经济理论认为,经济变量之间在水平上存在着均衡关系,可以 ...

  5. matlab泊松回归程序,R - 泊松回归( Poisson Regression)

    R - 泊松回归( Poisson Regression) 泊松回归涉及回归模型,其中响应变量是计数而不是分数的形式. 例如,足球比赛系列中的出生人数或获胜次数. 此外,响应变量的值遵循泊松分布. 泊 ...

  6. R语言用泊松Poisson回归、GAM样条曲线模型预测骑自行车者的数量

    全文链接:http://tecdat.cn/?p=18550 我根据泊松Poisson回归.GAM样条曲线模型对一个十字路口的骑自行车者的数量进行预测(点击文末"阅读原文"获取完整 ...

  7. R语言泊松回归(poisson)模型案例:基于robust包的Breslow癫痫数据集

    R语言泊松回归(poisson)模型案例:基于robust包的Breslow癫痫数据集 目录 R语言泊松回归(poisson)模型案例:基于robust包的Breslow癫痫数据集 #数据加载

  8. R语言惩罚逻辑回归、线性判别分析LDA、广义加性模型GAM、多元自适应回归样条MARS、KNN、二次判别分析QDA、决策树、随机森林、支持向量机SVM分类优质劣质葡萄酒十折交叉验证和ROC可视化

    最近我们被客户要求撰写关于葡萄酒的研究报告,包括一些图形和统计输出. 介绍 数据包含有关葡萄牙"Vinho Verde"葡萄酒的信息.该数据集有1599个观测值和12个变量,分别是 ...

  9. R语言主成分回归（PCR）、 多元线性回归特征降维分析光谱数据和汽车油耗、性能数据...

    原文链接:http://tecdat.cn/?p=24152 什么是PCR?(PCR = PCA + MLR)(点击文末"阅读原文"获取完整代码数据). • PCR是处理许多 x ...

  10. R语言构建logistic回归模型：构建模型公式、拟合logistic回归模型、模型评估，通过混淆矩阵计算precision、enrichment、recall指标

    R语言构建logistic回归模型:构建模型公式.拟合logistic回归模型.模型评估,通过混淆矩阵计算precision.enrichment.recall指标 目录

最新文章

  1. 12.5寸新i5商务本 联想X220i报10888元
  2. Apache 2.2 虚拟主机配置（本人推荐的）
  3. Delphi中高级DLL的编写和调用
  4. struts2找不到action_第一次用上Struts2框架做Web开发的体验……
  5. duilib入门简明教程 -- 部分bug (11) （转）
  6. 第十五期:真相了，中台到底“出路”还是“末路”?
  7. Java中对List集合排序的两种方法
  8. Win32 Thread Information Block
  9. plsql developer 无法登录Oracle
  10. 柴静：我只是讨厌屈服
  11. 每一代内存的读写速度
  12. python读书心得体会范文_读书心得体会范文5篇
  13. 玩冒险岛java卸载_冒险岛(经典版)卸载数据包方法
  14. 计算机白板培训心得,电子白板培训心得体会
  15. LeetCode第874题解析
  16. 《工程学导论》读书笔记第二章工程与科学
  17. ENVI 5.3 + 哨兵2号（Sentinel-2）L2A提取健康水体和不健康水体
  18. CF39C Moon Craters
  19. 统一身份认证和授权--微服务架构
  20. php开源广告系统,广告操作指南

热门文章

  1. 互联网的长在线、心跳和断线重连
  2. 应用程序的可视化图形引擎库控件Vectordraw Developer Framework
  3. mongodb安装指南 及使用
  4. (day9)357. 计算各个位数不同的数字个数
  5. solr4.2增量索引之同步(修改，删除，新增)
  6. Vuex的API文档
  7. 44个基于SaaS的商业智能解决方案
  8. Centos 6.5安装python3.5.1
  9. 编写安装配置DNS服务脚本
  10. android软键盘挡住输入框问题解决方法