【问题标题】:Binomial Test in SAS and R - Different ResultsSAS 和 R 中的二项式检验 - 不同的结果
【发布时间】:2019-01-01 09:44:24
【问题描述】:

我需要将二项式检验从 R 复制到 SAS,但我得到了不同的结果(或者我可能误解了 SAS 结果)。

为了简单地解释我的问题,我将使用来自这个维基百科的数据example,因为它提供了最终的解决方案;

假设您要计算在 235 个 6 面的公平骰子样本中获得 51 个或更多 6 的概率,因此每次试验中掷出 6 的真实概率为 1/6。 最终的解应该是0.02654

在 R 中,代码如下:

binom.test(51,235,(1/6),alternative = "greater")

得到的结果是:

精确二项式检验

数据:51 和 235 成功次数 = 51,试验次数 = 235,
p 值 = 0.02654
备择假设:成功的真实概率大于 0.1666667
95% 置信区间:
0.1735253 1.0000000
样本估计:成功的概率
0.2170213

在 SAS 中,等效项应为:

DATA DICEROLL;
ROLL=51;
FREQQ=235;
PROB=1/6;
RUN;

data _null_;
set diceroll;
call symput("probability",prob);
run;

PROC FREQ DATA=DiceRoll ;
    TABLES FREQQ / BINOMIAL (P=&probability.) ALPHA=0.05;
    EXACT  BINOMIAL ;
    WEIGHT ROLL ;
RUN;

但是THIS是我得到的结果(其中没有p值=0.02654)

我尝试了几种方法来协调我的结果(尝试了 R 中的所有三种替代方法,尝试在 sas 中反转 ROLL 和 FREQQ,因为我不确定)但我仍然没有找到解决方案。 binom.test 和 proc freq + BINOMIAL 是否至少执行相同的测试? 我是否误解了 SAS 输出?

提前感谢您的宝贵帮助!

=============================== 更新================ =============

我尝试了 reeza 和 BEMR 提出的两种方法,我觉得我已经接近解决方案了! @BEMR: 正如我在评论中所写和解释的那样,如果我的变量是二分的,我应该如何调整 %r(1,6) ?您的代码适用于 6 面骰子的示例,但在我的实际情况下,我的成功变量假定值介于 0 和 1 之间,所以我不确定我必须做什么(如果我没有在开始)

@REEZA:您的解决方案似乎有效,但我不得不删除 /2;我猜您的第一个解决方案将 p 值计算为双面测试而不是单面测试。 无论如何,结果都很好,但是当成功数为 0 或接近 0 (1,2,3) 时,SAS 和 R 之间存在巨大差异。您知道任何解决方法吗?或者更好的是,假设测试在这两种情况下都不可靠是否安全? 以下图片是我用reeza方法的结果,谢谢大家的宝贵合作!

【问题讨论】:

  • 看起来它认为您观察到了 100%,而不是 21.7% - 您如何指定成功次数?
  • 真的需要 SAS 中的 14 行来复制 R 中单行的功能吗?无论如何,n = 235 是一个相当大的样本量。您可以进行快速正态近似并独立验证 R 结果是否有意义,但任何表明 p 值
  • 不,它没有。在 R 中,您进行了理论计算,而在 SAS 中,您进行了模拟类型计算。而是在数据步骤中酌情使用 PPROB 或 QPROB 来获取概率。
  • 抱歉,您想要 CDF 或 PDF,但我的大脑现在有点炸了。我会发布我认为接近的内容,但您应该验证它。

标签: r sas binomial-theorem


【解决方案1】:

您显然不需要以这种方式设置变量,但这更像是一对一的类型比较。 SAS 没有能力进行我在函数中看到的单方面测试,但我没有阅读太多内容或尝试确定它是否正确。但这是您应该在 SAS 中使用的方法类型来获得相似的数字,而不是 PROC FREQ。

    data demo;
nSuccesses=51;
prob_success=1/6;
nTrials = 235;

y=(1-cdf('BINOM', nsuccesses, prob_success, ntrials))/2;
run;

proc print data=demo;
run;

http://documentation.sas.com/?docsetId=lefunctionsref&docsetTarget=p1cxa81efqtsszn12ueyitll9esw.htm&docsetVersion=9.4&locale=ja#p03dt2kdzjjucxn198ytlpnrf1r4

【讨论】:

  • 通过编辑您的问题发布您拥有的代码,我会再看看它。
  • 对不起,今天我也有点炸了!我的评论提到了我尝试使用的另一种方法,我使用它获得了置信区间但没有 p 值 (1-betainv(.025,(&trials.-&nsuccess.),&nsuccess.+1)。您的方法似乎有效,但效果不大问题。我会尽快更新我的问题
【解决方案2】:

如果您想比较 binom.test 和 proc freq + BINOMIAL,您可以在 SAS 中使用模拟。以下代码提供了一个示例:

当骰子掷出 235 时,结果可能是 1,...,6。

*Create df: random roll;
*macro: random int between min and max;;
  %macro r(min,max);
(&min + floor((1+&max-&min)*rand("uniform"))) 
   %mend;
  data df;
  f = 0;
  do i = 1 to 235; *number of trials;
    x = %r(1,6); *call macro %r() to generate random number between 1,...,6; 
if x = 6 then f = f + 1; *if the random number = 6, add freq from the previous;
relative = f/i; *relative freq;
 output;
end;
run;

*plot relative freq, reference line (1/6), probability of rolling 6;
symbol v=dot c=red;
proc gplot data=df;
plot relative * i/overlay vref=0.16666667 href=500 lh=3;
run;
quit;

以下是此处的示例:http://www.stat.purdue.edu/~lfindsen/stat503/Lab2.pdf

*exact binomial using proc freq and simulated data; 
*test if simulation is different from the hypothized 1/6;
proc freq data = df;
tables x / binomial (level=6 p=.166667); 
exact binomial;
run;

当 51 例是 235 例中的 6 例时。

*Create df2: assign approx 51 cases of 235 a roll of 6;
 data df2;
 do i = 1 to 235; *number of trials;
x = %r(1,5); 
 output;
 end;
 run;
 data df2;
  set df2;
    if i <= 51 then x = 6; *assign six to rows 1 to 51; 
  run;

  *exact binomial using proc freq and simulated data; 
  *test if simulation is different from the hypothized 1/6;
  proc freq data = df2;
tables x / binomial (level=6 p=.166667); 
exact binomial;
  run;

精确的二项式单边 p 值 = 0.0265

===============================更新============= ===============

对于二进制变量 [0=2184,1=72] 而不是使用宏,您可以执行以下操作:

    data df3;
    input success n;
    datalines;
    0 2184
    1 72
    ;

    proc freq data=df3;
    weight n; *number of obs for [0,1];
    tables success / binomial (level=2 p=0.509); 
    run;

【讨论】:

  • 我需要做的是检查估计的概率是否低估了观察到的成功率。假设你有:Success=72;试验=2256;成功率=3.19%;概率估计=5.09%。实验只能假设值 0(失败)或 1(成功)。 R 中的零假设是成功的真实概率小于 0.509,因此对于 p 值 > alpha,测试通过。 R 代码将是:> binom.test(72,2256,0.509,alternative = "greater")。我需要的结果 -> p-value = 1。如果我的变量是二分的,我应该如何调整 %r(1,6)?
  • 如果我得到了正确的代码逻辑,它应该是 %r(0,0) 然后如果 i
  • @BEMR 为什么在测试双面时,SAS 的双面 p 值为 0.0531?使用 R,你会得到 0.04375 的 p 值?
猜你喜欢
  • 1970-01-01
  • 2021-01-21
  • 2017-11-09
  • 2015-09-24
  • 2017-07-05
  • 2021-10-22
  • 1970-01-01
  • 2019-06-02
  • 2015-04-09
相关资源
最近更新 更多