首页 > 其他分享 >R语言进行数值模拟:模拟泊松回归模型的数据

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

时间:2024-01-22 21:55:34浏览次数:18  
标签:泊松 系数 方差 xb 模型 数值 回归 模拟

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

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

 

模拟回归模型的数据

验证回归模型的首选方法是模拟来自它们的数据,并查看模拟数据是否捕获原始数据的相关特征。感兴趣的基本特征是平均值。我喜欢这种方法,因为它可以扩展到广义线性模型(logistic,Poisson,gamma,...)和其他回归模型,比如t -regression。

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

我将用一个例子来证明我的意思。

示范

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

   
n <- 50
set.seed(18050518) 
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

与忽略系数不确定性时相比,均值和方差更接近原始结果的均值和方差。与考虑回归系数的不确定性时相比,这种方法总是会导致方差较小。它要快得多,并且需要零编程来实现,但我不习惯忽略回归系数的不确定性,使模型看起来比它更充分。

 

非常感谢您阅读本文,有任何问题请在下面留言!

最受欢迎的见解

1.R语言多元Logistic逻辑回归 应用案例

2.面板平滑转移回归(PSTR)分析案例实现

3.matlab中的偏最小二乘回归(PLSR)和主成分回归(PCR)

4.R语言泊松Poisson回归模型分析案例

5.R语言回归中的Hosmer-Lemeshow拟合优度检验

6.r语言中对LASSO回归,Ridge岭回归和Elastic Net模型实现

7.在R语言中实现Logistic逻辑回归

8.python用线性回归预测股票价格

9.R语言如何在生存分析与Cox回归中计算IDI,NRI指标

 

标签:泊松,系数,方差,xb,模型,数值,回归,模拟
From: https://www.cnblogs.com/tecdat/p/17981170

相关文章

  • 2024 蓝桥杯模拟赛 1 (div1) 题解
    A.把字符串小写转换成大写即可#include<bits/stdc++.h>usingnamespacestd;voidsolve(){strings;cin>>s;for(inti=0;i<s.size();i++){if(s[i]>='a'&&s[i]<='z'){s[i]=(char)(s[i]-'a......
  • 路由器固件模拟环境搭建
    路由器固件模拟环境搭建binwalk安装参考参考链接https://xz.aliyun.com/t/5697?time__1311=n4%2BxnD07Dti%3D0%3DDk8GCDlhjm5fcQQeiKN4D&alichlgref=https%3A%2F%2Fwww.google.com%2Fhttp://zeroisone.cc/2018/03/20/固件模拟调试环境搭建/但是他们都有一个问题,在按他们的步......
  • 【快速阅读三】使用泊松融合实现单幅图的无缝拼贴及消除两幅图片直接的拼接缝隙。
    泊松融合还可以创建一些很有意思的图片,比如一张图片任意规格平铺,使用泊松融合后,平铺的边界处过渡的很自然,另外,对于两张图片,由于局部亮度等等的影响,导致拼接在一起时不自然,也可以使用泊松融合予以解决。在【快速阅读二】从OpenCv的代码中扣取泊松融合算......
  • 【快速阅读二】从OpenCv的代码中扣取泊松融合算子(Poisson Image Editing)并稍作优化
    泊松融合是一种非常不错多图融合算法,在OpenCv的相关版本中也包含了该算子模块,作者尝试着从OpenCv的大仓库中扣取出该算子的全部代码,并分享了一些在扣取代码中的心得和收获。泊松融合我自己写的第一版程序大概是2016年在某个小房间里折腾出来的,当时是用......
  • 2024.1.21模拟赛 C题解
    简要题意略思路首先有一个\(O(nk)\)的暴力dp,30pts我们可以发扬人类智慧,构造势能函数\(U_x=\sum_{未选择的点i}dis(i,x)+h_i\),当前在\(x\)点定义\(f_i\)表示走到\(i\)点时势能函数的最小值,\(s_i\)表示\(i\)到起点的距离容易发现只会跨过起点进行转移,于是\(f_i=f_j+2\tim......
  • 2024.1.21模拟赛 B题解
    题目大意略思路首先有一个50pts的网络流暴力考虑按照\(dp\)值分层,发现在同一层内,随着\(i\)递增,\(a_i\)递减由此可以进一步推出每一个点连接的出边,是下一层的一个区间,并且区间是单调的于是可以线段树优化建边,拿到60pts接着考虑模拟网络流,发现如果每次都选择第一条出边的话,就......
  • 初三选拔模拟赛题解
    目录T2T3T4T5T6T7给个正常的题解以正视听.不过说好的普及难度呢?如果有问题请指出.T2注意到答案一定可以取到最小区间的长度\(len\),一种方案是按\(0\dotslen-1\)循环填.T3大致有两种做法:维护每个手指的次数\(c_i\)和占用的键数\(t_i\),按\(\frac{c_i}{t_i+1}\)......
  • 1.21 && 第二场模拟赛记
    写在前面:非常好模拟赛,9道题,3道不用写,三道原题,两道原题,一道东方题。根据等量代换可得有5道原题。t2原题CF740C赛时理解错题意了,具体咋想的我也忘了,但是我的构造方法是每个区间从0开始构造,如果不在区间内则任意输出。但是正解是找到最小区间然后按区间循环输出即可。......
  • 模拟扑克牌初始化操作
    1publicclassPoker{2/*模拟扑克牌初始化操作3*点数:"3","4","5","6","7","8","9","10","J","Q","K","A","2","大......
  • 1.21闲话暨模拟赛题解
    未卜先知推歌:二重变革/洛天依,言和byDELA上午写了模拟赛,下午不给我发代码不让我改题不让我看题面,smjb模拟赛一共9道题,4道原题(2道原题,2道"原"题),抛去3道不写的一共6道题T1「尘世闲游」(原神题)没让写T2「一心净土」(原神题&&原题「CF740C」)题解我这里一共找到......