R语言进行数值模拟:模拟泊松回归模型的数据

验证回归模型的首选方法是模拟来自它们的数据,并查看模拟数据是否捕获原始数据的相关特征。

由Kaizong Ye,Coin Ge撰写

感兴趣的基本特征是平均值。我喜欢这种方法,因为它可以扩展到广义线性模型(logistic,Poisson,gamma,…)和其他回归模型,比如t -regression

模拟回归模型的数据

标准回归模型假设存在将预测变量与结果相关联的真实/固定参数。但是,当我们执行回归时,我们只估计这些参数。因此,回归软件返回表示系数不确定性的标准误差。

×

泊松回归很显然是y符合泊松分布或者假设y符合泊松分布(到底是不是真的泊松分布并不重要),对应“广义线性模型”中的泊松分布模型。如果“链接函数”是g(u)=u(这种链接函数用的非常少),则泊松回归并不属于“对数线性模型”,只有g(u)=ln(u)时才属于“对数线性模型”。

“对数线性模型” 包含了很多模型,其数学模型常见的等式是:ln(u)=a+X’B,a是一个常量,X是自变量向量,B是回归系数向量。注意这个等式的右边是没有误差项的,为什么不需要误差项呢?因为u是一个确定值,不是随机变量。u可以是某个分布(比如泊松分布,负二项分布之类)的期望,或者一个其他的确定函数(比如生存分析中的风险率)。至于模型是否可以有误差项,我也不是很清楚。毕竟“对数线性模型”这个模型不是特指,是一个泛指。

经典线性回归模型是 y=a+X’B+w ,y是随机变量,w是误差项,是随机变量。有时候为了解决异方差、y的分布是对数正态、y是右偏的或者y的取值太分散,会对y取对数,模型变成 ln(y)=a+X’B+w。这个模型到底算不算“对数线性模型”呢,应该也算吧。但是个人感觉把它作为经典线性模型看待可能更好点,因为这个对数线性模型使用了经典线性模型的假定和求解方法。

对数线性模型使用最多的可能是列联表分析,列联表的每个格子数目可以看成符合泊松分布,或者所有格子数符合多项分布。列联表某个格子的期望是u=aAiBj,这个乘积取对数后就成了“对数线性模型”。

生存分析里面的“cox回归” 模型是:h(t|X)=h0(t)g(X),假设g(X)是指数函数,取对数后模型成了“对数线性模型”。

Logistic 模型算不算“对数线性模型” 呢?好像也可以算,不过它已经有了自己的专有名字,所以一般不叫它们是“对数线性模型” 。

社会学,生物统计学和医学统计学里面说“对数线性模型”通常指的是分析列联表的“对数线性模型”,比如SPSS软件里面的“对数线性模型”模块。Stata 软件可视化界面里面没有“对数线性模型”这个模块(不过可以使用命令行做列联表的对数线性模型分析),可能是stata 在计量经济学用的比较多。


我将用一个例子来证明。

示范

我将使用泊松回归来证明这一点。我模拟了两个预测变量,使用50的小样本。

n <- 50
set.seed(18050518) 


课程

R语言数据分析挖掘必知必会

从数据获取和清理开始,有目的的进行探索性分析与可视化。让数据从生涩的资料,摇身成为有温度的故事。

立即参加

xc的系数为0.5 ,xb的系数为1 。我对预测进行取幂,并使用该rpois()函数生成泊松分布结果。

#指数预测并传递给rpois()

summary(dat)

       xc                  xb             y       
 Min.   :-2.903259   Min.   :0.00   Min.   :0.00  
 1st Qu.:-0.648742   1st Qu.:0.00   1st Qu.:1.00  
 Median :-0.011887   Median :0.00   Median :2.00  
 Mean   : 0.006109   Mean   :0.38   Mean   :2.02  
 3rd Qu.: 0.808587   3rd Qu.:1.00   3rd Qu.:3.00  
 Max.   : 2.513353   Max.   :1.00   Max.   :7.00  

接下来是运行模型。

 
Call:
glm(formula = y ~ xc + xb, family = poisson, data = dat)

Deviance Residuals:
    Min       1Q   Median       3Q      Max  
-1.9065  -0.9850  -0.1355   0.5616   2.4264  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.20839    0.15826   1.317    0.188    
xc           0.46166    0.09284   4.973 6.61e-07 ***
xb           0.80954    0.20045   4.039 5.38e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 91.087  on 49  degrees of freedom
Residual deviance: 52.552  on 47  degrees of freedom
AIC: 161.84

Number of Fisher Scoring iterations: 5

估计的系数与人口模型相距不太远,.21代表截距而不是0,.46而不是.5,而0.81而不是1。

接下来模拟模型中的数据,我想要10,000个模拟数据集,为了捕捉回归系数的不确定性,我假设系数来自多元正态分布,估计系数作为均值,回归系数的方差 – 协方差矩阵作为多元正态分布的方差 – 协方差矩阵。

coefs <- mvrnorm(n = 10000, mu = coefficients(fit.p), Sigma = vcov(fit.p))

检查模拟系数与原始系数的匹配程度。

coefficients(fit.p)

(Intercept)          xc          xb
  0.2083933   0.4616605   0.8095403

colMeans(coefs) #模拟系数的平均值

(Intercept)          xc          xb
  0.2088947   0.4624729   0.8094507

标准错误:

sqrt(diag(vcov(fit.p)))

(Intercept)          xc          xb
 0.15825667  0.09284108  0.20044809

apply(coefs, 2, sd) #模拟系数的标准偏差

(Intercept)          xc          xb
 0.16002806  0.09219235  0.20034148

下一步是模拟模型中的数据。我们通过将模拟系数的每一行乘以原始预测变量来实现。然后我们传递预测:

#每种情况一行,每组模拟系数一行

#模型矩阵与系数的乘积,取幂,
#然后用于模拟泊松分布的结果
for (i in 1:nrow(coefs)) {
  sim.dat[, i] <- rpois(n, exp(fit.p.mat %*% coefs[i ]))
}
rm(i, fit.p.mat) # 

现在一个是完成模拟,将模拟数据集与原始数据集至少比较结果的均值和方差:

c(mean(dat$y), var(dat$y))#原始结果的均值和方差

[1] 2.020000 3.366939

c(mean(colMeans(sim.dat)), mean(apply(sim.dat, 2, var)))#10,000个模拟结果的平均值和变量的平均值

[1] 2.050724 4.167751

模拟结果的平均值略高于原始数据,平均方差更高。平均而言,可以预期方差比平均值更偏离目标。方差也将与一些极高的值正偏差,同时,它的界限为零,因此中位数可能更好地反映了数据的中心:

 
[1] 3.907143

中位数方差更接近原始结果的方差。

这是模拟均值和方差的分布:


 绘制10,000个模拟数据集中的一些数据集的直方图并将其与原始结果的直方图进行比较也是有用的。还可以测试原始数据和模拟数据集中xb = 1和xb = 0之间结果的平均差异。 


回到基础R,它具有simulate()执行相同操作的功能:

sim.default <- simulate(fit.p, 10000)

此代码相当于:

sim.default <- replicate(10000, rpois(n, fitted(fit.p)))

fitted(fit.p)是响应尺度的预测,或指数线性预测器,因为这是泊松回归。因此,我们将使用模型中的单组预测值来重复创建模拟结果。

 c(mean(colMeans(sim.default)), mean(apply(sim.default, 2, var)),
 
[1] 2.020036 3.931580 3.810612

与忽略系数不确定性时相比,均值和方差更接近原始结果的均值和方差。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498

关注有关新文章的微信公众号


永远不要错过任何见解。当新文章发表时,我们会通过微信公众号向您推送。

技术干货

最新洞察

This will close in 0 seconds