【統計】時系列分析（ARIMA、SARIMA、SARIMAXをSASで実行する）

SAS 統計

ARIMAプロシジャのメモ。
基本はRを使っているので初歩。

【目次】


以前の記事
【統計】時系列分析（AR、MA、ARMA、ARIMA、SARIMA、ARIMAXモデル） - こちにぃるの日記
【統計】時系列分析（ARIMA、SARIMA、SARIMAXをRで実行する） - こちにぃるの日記


Auto ARIMA

identify ステートメントでAR, MA過程の複数パターンを評価できる。
var=target で目的変数を指定。カッコ内に階差を指定。
nlag で季節サイクルを指定。
p=(0:10) q=(0:10) でAR過程を０～10の各組合せで推定。
stationarity=(adf) で定常性を検定。adf で拡張ディッキー–フラー検定。
crosscorr=(variable_name) で外生変数（回帰成分）との相関を評価。
esacf minic scan あたりで各組合せを評価。

proc arima data=adts ;
    identify var=target(1,12) nlag=12
             p=(0:10) q=(0:10)
             stationarity=(adf)
             crosscorr=(variable_name) outcov=_outcov
             esacf minic scan ;
run;
quit;

ARIMA(p,0,0)：AR(p)モデル

estimate ステートメントでAR過程とMA過程を指定。

proc arima data=adts ;
    identify var=target;
    estimate p=2; 
run;
quit;

ARIMA(0,0,q)：MA(q)モデル

proc arima data=adts ;
    identify var=target;
    estimate q=2; 
run;
quit;

ARIMA(p,0,q)：ARMA(p,q)モデル

proc arima data=adts ;
    identify var=target;
    estimate p=2 q=2; 
run;
quit;

ARIMA(p,d,q)：ARIMA(p,d,q)モデル

proc arima data=adts ;
    identify var=target(1);
    estimate p=2 q=2; 
run;
quit;

ARIMAXモデル

identify ステートメントの crosscorr と estimate ステートメントの input で外生変数を指定。

proc arima data=adts ;
    identify var=target(1) crosscorr=(variable_name);
    estimate p=2 q=2 input=(variable_name);
run;
quit;

SARIMAモデル

nlag でサイクルを指定。

proc arima data=adts ;
    identify var=target(1,12) nlag=12;
    estimate p=(1 2)(1 2) q=(1 2)(1 2); 
run;
quit;

SARIMAXモデル

proc arima data=adts ;
    identify var=target(1,12) nlag=12 crosscorr=(variable_name);
    estimate p=(1 2)(1 2) q=(1 2)(1 2) input=(variable_name);
run;
quit;

将来予測

forecast ステートメントで将来予測。
lead で何期まで予測するか指定。
interval でサイクルの単位。
id で inputデータの必要な変数を outputデータに保持させる。

proc arima data=adts ;
    identify var=target(1,12) nlag=12 crosscorr=(variable_name);
    estimate p=(1 2)(1 2) q=(1 2)(1 2) input=(variable_name);
    forecast lead=12 interval=month id=id_variable_name out=_pred ;
run;
quit;
  
proc sgplot data=_pred;
   band Upper=u95 Lower=l95 x=id_variable_name / LegendLabel="95%CI";
   scatter x=id_variable_name y=target;
   series  x=id_variable_name y=forecast ;
run;


参考

https://support.sas.com/documentation/onlinedoc/ets/132/arima.pdf

【統計】割合の信頼区間（WaldとClopper-pearson）

R、R Studio SAS 統計

割合の信頼区間について

【目次】


計算式

Wald (正規近似を用いた方法)
\quad \widehat{p}\pm z\sqrt{\dfrac{\widehat{p}\left( 1-\widehat{p}\right) }{n}}

二項分布 B (n, p) は n が十分大きいとき、平均 np , 分散 np(1-p) の正規分布に近づく（ラプラスの定理）。 一般的に、np>5 かつ n(1-p)>5 の場合を n が十分に大きいと考えるようである。

n が十分に大きくない場合は近似が悪い。また、０や１の両端も近似が悪いので不均衡データ（Imbalanced Data）に注意が必要となる。


Clopper-pearson (F分布を用いた方法)
\quad Lower = \Biggl( 1+ \dfrac{n-x+1}{(x)F\bigl( 1-\frac{\alpha}{2}; \ 2(x) , \ 2(n-x+1) \bigr) } \Biggr)^{-1}
\quad Upper = \Biggl( 1+ \dfrac{n-x}{(x+1)F\bigl( \frac{\alpha}{2}; \ 2(x+1) , \ 2(n-x) \bigr) } \Biggr)^{-1}

二項分布（離散確率分布）をベータ分布（連続値型確率分布）に拡張し、F分布で表したものを利用した信頼区間。

計算例

Wald
=D4-NORM.S.INV(1-G4/2)*SQRT(D4*(1-D4)/C4)
=D4+NORM.S.INV(1-G4/2)*SQRT(D4*(1-D4)/C4)

Clopper-pearson
=(1+(C4-B4+1)/(B4*F.INV.RT(1-G5*0.5, 2*B4, 2*(C4-B4+1))))^-1
=(1+(C4-B4)/((B4+1)*F.INV.RT(G5*0.5, 2*(B4+1), 2*(C4-B4))))^-1


コード

library(DescTools)
DescTools::BinomCI(x, n, conf.level=.95,method=c("wald"))
DescTools::BinomCI(x, n, conf.level=.95,method=c("clopper-pearson"))



proc freq data=ADS;
    tables YN / binomial;
    weight CNT / zeros;
run;


被覆確率 (Coverage Probability)

被覆確率（信頼区間が真の値を含む確率）を調べてみる。
p = ０ ～ １を 0.01 刻みで、二項分布 B(n, p) に従う乱数を100個ずつ生成（PCのスペックがないのでこのくらいで）。
n = 10, 100, 1000 の３種類を試してみる。

被覆確率 (Coverage Probability) をグラフ化し、95％のところに横線を引いている。
n が小さいとき、Waldは被覆確率95％を大きく下回る。両端の不均衡データにあたる部分でその傾向が大きく見られる。
Clopper-pearsonは n が小さくても95%から大きく下回ることはなさそう。



参考

https://www.lexjansen.com/phuse/2013/sp/SP05.pdf
https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups14-d-05.pdf
統計学入門−第3章

【統計】割合の差の信頼区間（リスク差の信頼区間）

R、R Studio SAS 統計

SASとRでの使い方メモ。

　 R
DescTools::BinomDiffCI
method = TYPE		 SAS
PROC FREQ
RISKDIFF (CL = TYPE)
Agresti-Caffo  "ac" AC
Brown, Li's Jeffreys "blj" -
Farrington and Manning - FM
Haldane "hal" -
Hauck-Anderson "ha" HA
Jeffreys-Perks "jp" -
Mee  "mee" MN (CORRECT= NO)
Miettinen and Nurminen  "mn" MN
Newcombe (Corrected) "scorecc" NEWCOMBE (CORRECT)
| SCORE (CORRECT)
| WILSON (CORRECT)
Newcombe  "score" NEWCOMBE
| SCORE
| WILSON
Wald (Corrected) "waldcc" WALD (CORRECT)
Wald  "wald" WALD


Rでの実行

c00 <- 50; c01 <- 75; N0 <- c00+c01;
c10 <- 40; c11 <- 90; N1 <- c10+c11;
 
library(DescTools)
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="ac")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="ha")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="mee")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="mn")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="score")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="scorecc")
100 * BinomDiffCI(c00, N0, c10, N1, conf.level=.95, method="wald")


SASでの実行

data ADS;
    length FACTOR YN CNT 8.;
    input FACTOR YN CNT @@;
    cards;
0 0 50  0 1 75  1 0 40  1 1 90
;
run;
  
proc freq data=ADS;
    table FACTOR * YN / riskdiff(cl=(wald ha fm ac mn newcombe));
    table FACTOR * YN / riskdiff(cl=(wald(correct) mn(correct=no) newcombe(correct)));
    weight CNT / zeros;
run;


参考

https://cran.r-project.org/web/packages/DescTools/DescTools.pdf
https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups14-d-08.pdf
https://www.sas.com/content/dam/SAS/ja_jp/doc/event/sas-user-groups/usergroups2015-b-01.pdf
https://www.lexjansen.com/wuss/2016/127_Final_Paper_PDF.pdf
https://www.seikei.ac.jp/university/rikou/pdf/V0410202.pdf
https://www.ms.uky.edu/~mai/sta635/FagerlandLydersenLaake2011---RecommendedCIsForTwoIndependent....pdf
https://onlinelibrary.wiley.com/doi/epdf/10.1002/pst.2162

【統計】ロジスティック回帰のFisher Scoringによる最適化

R、R Studio SAS 統計 機械学習


【目次】


以下でロジスティック回帰に触れていたがFisher Scoringにて再び。
【統計】ロジスティック回帰分析 - こちにぃるの日記
【統計】ロジスティック回帰分析（ホスマー・レメショウ検定、ROC、AUC等） - こちにぃるの日記
【統計】一般化線形モデル（GLM） - こちにぃるの日記
【Python】Kerasでロジスティック回帰 - こちにぃるの日記
【Kaggle】タイタニックの分析をやってみる - こちにぃるの日記


ロジスティック回帰の特徴

ロジスティック回帰は一般化線形モデルの一つ。

目的変数 二値の離散値
分布族 二項分布
リンク関数 ロジットリンク関数
逆リンク関数 ロジスティック関数（シグモイド関数）

・誤差構造が二項分布に従うと仮定。
・ロジット関数で実測値 (０～１) をロジット値 (-∞～∞) に変換し、線形モデルを作成。

\qquad logit(θ_{i})=ln (\dfrac{θ_{i}}{1-θ_{i}})=X_{i}' \beta

・ロジスティック関数で線形予測子から予測値（０～１）を得る。

\qquad θ_{i} = logistic( X_{i}' \beta )=\dfrac{1}{1+exp( -X_{i}' \beta )}=\dfrac{exp( X_{i}' \beta )}{1+exp( X_{i}' \beta )}

・二値データに対する通常の回帰モデルでは、予測値が０～１の範囲を超えてしまい解釈が難しいが、 ロジスティック回帰の予測値は０～１の範囲に収まり、解釈しやすい。

（青線：ロジスティック回帰、緑線：単回帰分析）


ロジスティック回帰の尤度・誤差関数

尤度

\quad L \left( \beta \right) =\Pi \{ θ_{i}^{\ y_{i}} \times \left( 1-θ_{i} \right)^{ \left( 1 - y_{i} \right) } \}

対数尤度

\quad ln\ L\left( \beta\right)=\sum \{ y_{i} \cdot ln\left(θ_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-θ_{i} \right) \}

交差エントロピー誤差（分類問題の損失関数）

\quad E\left(\beta\right)=-\sum \{ y_{i} \cdot ln\left(θ_{i}\right) + \left(1-y_{i} \right) \cdot ln \left( 1-θ_{i} \right) \} =0


下図の通り、同一式で Y=１ or ０ の実測値と予測値の差（誤差）が求まる。


ロジスティック回帰のFisher Scoringによる最適化

勾配

\quad \nabla E( \beta )=\dfrac{\partial E( \beta) }{\partial θ}\dfrac{\partial θ}{\partial \widehat{y}}\dfrac{\partial \widehat{y}}{\partial \beta }

\qquad \qquad =-\sum \{ y_{i} \dfrac{1}{θ_{i}} θ_{i}( 1-θ_{i}) x_{i} + (1-y_{i}) \dfrac{1}{1-θ_{i}} ( -θ_{i}( 1-θ_{i} ) ) x_{i}  \}

\qquad \qquad =-\sum \{ y_{i} ( 1-θ_{i}) + (1-y_{i}) ( -θ_{i} )\} x_{i}

\qquad \qquad =-\sum \{ y_{i} - y_{i} θ_{i} - θ_{i} + y_{i} θ_{i}\} x_{i}

\qquad \qquad =-\sum \{ y_{i} - θ_{i} \} x_{i}

\qquad \qquad =\sum ( θ_{i} - y_{i}) x_{i}=0

\quad \nabla E( \beta )= X' (θ-Y)


Fisher情報行列

\quad I=-E(\nabla^{2} E( \beta ))

\quad \nabla^{2} E( \beta ) =-\sum ^{n}_{i=1}\{x_{ik}' \ \theta _{i}(1-\theta _{i}) \ x_{ik})\}

\quad \nabla^{2} E( \beta ) =-X'WX


更新式（Fisher Scoring）
Hessian行列を期待値に置き換えたFisher Scoring法。

\quad \beta^{(t+1)}=\beta^{(t)}-(I^{(t)})^{-1}\nabla E( \beta )^{(t)}

\qquad\quad\ = \beta^{(t)} - (X'WX)^{-1} X'(θ-Y)


手計算

Kaggleチュートリアルのタイタニックデータを用いる。
簡単に以下のモデルとする。

目的変数 Survived（生存有無） ０= 死亡, １= 生存
説明変数 Pclass（チケットのクラス） １= 1st, ２= 2nd, ３= 3rd
※今回は簡単に ０= 1st or 2nd, １= 3rd
Sex（性別） ０= Female, １= Male

Rでの実行結果

Fisher Scoring iterations は 4回。

titanic$Event   <- 1 - titanic$Survived
titanic$Pclass3 <- if_else(titanic$Pclass=="3",1,0)
titanic$Sex2    <- if_else(titanic$Sex=="male",1,0)
  
Call:
glm(formula = Event ~ Pclass3 + Sex, family = binomial, data = titanic)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.1603  -0.5310   0.4516   0.8910   2.0143  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.8878     0.1776 -10.629   <2e-16 ***
Pclass3       1.5125     0.1783   8.483   <2e-16 ***
Sexmale       2.6067     0.1818  14.337   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1186.66  on 890  degrees of freedom
Residual deviance:  838.86  on 888  degrees of freedom
AIC: 844.86

Number of Fisher Scoring iterations: 4


手計算
初期値は経験ロジットによる回帰パラメータを与える（LINEST関数で推定）。

【更新１回目】

【更新２回目】

【更新３回目】

【更新４回目】

更新４回目で回帰パラメータの前後に差がなくなったため終了する。
推定パラメーターはRと同じ。

β exp(β)
X0 -1.88780 0.15140
X1 1.51251 4.53810
X2 2.60674 13.55479

チケットが1st/2ndクラスに比べて3rdの方が4.5倍オッズが高い。
女性に比べて男性の方が13.6倍オッズが高い。


参考

https://ocw.tsukuba.ac.jp/wp/wp-content/uploads/2019/10/5b43d44b7c8f6475622e63173d226184.pdf
newton
Estimating Logistic Regression Coefficents From Scratch (R version) - The Pleasure of Finding Things Out
Estimating Logistic Regression Coefficents From Scratch (Python version) - The Pleasure of Finding Things Out
http://homepage.stat.uiowa.edu/~luke/classes/STAT7400-2021/_book/optimization.html

【統計】ウィルコクソンの符号付き順位検定（Exact）

R、R Studio SAS 統計

Wilcoxonの符号付き順位検定のExact法（正確確率検定）のメモ。

【目次】


正規近似とT近似での検定は以下記事。
cochineal19.hatenablog.com


計算例

簡単なデータを準備。
下図は正規近似・T近似での符号付き順位検定の記事で使った表。
ここにSinged Rank列（符号付き順位、前後差がマイナスなら順位にマイナスを付ける）とSinged Rank列の合計を追加（紫色のセル）。


観測値（前後差の絶対値）の順位（今回は１、2.5、2.5、４、５位）を固定した場合に、符号付き順位がとり得る全パターンとその順位和を計算する （ゼロ処理はWilcoxon法：前後差ゼロは含めない、タイは平均順位）。 下図の25=32パターンが考えられる。


上図の全パターンを順位和別に頻度集計し、生起確率を計算する。
今回、観測値の符号付順位和は７なので、順位和=７とそれより珍しいパターンの確率を足し合わせる。
順位和=７～15が該当し、この確率の和が片側確率、その２倍が両側確率となる。


Rで実行。
coinパッケージの wilcoxsign_test を使用。
zero.method="Wilcoxon" でWilcoxon法を指定。
distribution="exact" で正確確率を指定。

> ads0 <- data.frame(
+   before=c(10,13,13,13,20,21,21)
+   ,after=c(19,18,18,15,14,21,21)
+ )
> coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="exact")

    Exact Wilcoxon Signed-Rank Test

data:  y by x (pos, neg) 
     stratified by block
Z = -0.9482, p-value = 0.4375
alternative hypothesis: true mu is not equal to 0


プログラムコード

library(coin)
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="exact")



proc univariate data=ads;
    var diff;
    output out=out1 probs=probs;
    ods output TestsForLocation=out2;
run;

SASでは、n<20はExact（正確確率）、n>=20ならWilcoxon法・t近似のp値を返す。


参考

SAS Help Center
符号検定とWilcoxonの符号付き順位検定

【統計】Wilcoxonの順位和検定（Exact）

SAS R、R Studio 統計

Wilcoxonの順位和検定のExact法（正確確率検定）のメモ。

【目次】


正規近似とT近似での検定は以下記事。
cochineal19.hatenablog.com


計算例

簡単なデータを準備。
A群が３例で順位和=７、B群が４例で順位和=21。


観測値の順位（今回は１～７位）を固定した場合に、A群が配置されるパターンとその順位和を計算する。
A群３例の配置は下図の35パターンが考えられる（A群=１、B群=０）。


上図の全パターンを順位和別に頻度集計し、生起確率を計算する。
今回、A群の観測値の順位和は７なので、順位和=７とそれより珍しいパターンの確率を足し合わせる。
順位和=６と７が該当し、この確率の和が片側確率、その２倍が両側確率となる。

なお、正規近似での検定は、上記ヒストグラムに正規分布を当てはめるイメージ。


Rで実行。
coinパッケージの wilcox_test を使用。distribution="exact" で正確確率を指定。

> ads <- data.frame(
+   Group <- factor(c("A","A","B","A","B","B","B"))
+   ,AVAL <- seq(10:16))
> wilcox_test(AVAL ~ Group, data=ads, method="exact", distribution="exact")

    Exact Wilcoxon-Mann-Whitney Test

data:  AVAL by Group (A, B)
Z = -1.7678, p-value = 0.1143
alternative hypothesis: true mu is not equal to 0


SASではexactステートメントで指定する。

data ads;
  do AVAL = 10,11,13; Group=1; output; end;
  do AVAL = 12,14,15,16; Group=2; output; end;
run;
proc npar1way data=ads wilcoxon correct=no /*連続修正なし*/;
  class Group;
  var   AVAL;
  exact wilcoxon;
run;


プログラムコード

library(coin)
coin::wilcox_test(AVAL ~ Group, data=ads, method="exact", distribution="exact")



proc npar1way data=ads wilcoxon correct=no ;
  class TRT;
  var   AVAL;
  exact wilcoxon;
run;


参考

Wilcoxon Rank Sum Exact Test | Real Statistics Using Excel

【統計】ウィルコクソンの符号付き順位検定（Wilcoxon Signed Rank Test）

SAS R、R Studio

ウィルコクソンの符号付き順位検定（Wilcoxon Signed Rank Test）のメモ。
ノンパラメトリックな対応のある差の検定（２時点）。

【目次】

正規近似、T近似について。exactは記載なし。

帰無仮説、対立仮設

・帰無仮説 H_{0}: ２時点に差がない

・対立仮設 H_{1}: ２時点に差がある


計算式等

■ 計算の流れ
① 前後差（変化量）を算出する。
② 前後差の絶対値に平均順位を付ける。前後差ゼロの扱いは手法により異なる。
　 Wilcoxon法：ゼロを除いて平均順位を付ける。
　 Pratt法　　：ゼロを含めて平均順位を付ける。
③ 前後差が正の平均順位和 T^{+} を用いて検定統計量 S を算出し、p値を得る。


■ 期待値・分散
　< Wilcoxon法 >

\quad E(T)=\dfrac{n(n+1)}{4}

\quad V(T)=\dfrac{n(n+1)(2n+1) - 0.5 \sum t_{i}(t_{i}+1)(t_{i}-1)}{24}


　< Pratt法 >

\quad E(T)=\dfrac{n(n+1) -t_{0}( t_{0}+1)}{4}

\quad V(T)=\dfrac{n(n+1)(2n+1) - t_{0}(t_{0}+1)(2t_{0}+1) - 0.5 \sum t_{i}(t_{i}+1)(t_{i}-1)}{24}


　※ n=被験者数、t_{i}=前後差の絶対値のタイ数、t_{0}=前後差ゼロの数


■ 統計量

\quad S=T^{+} - E(T)

　※ T^{+}=前後差が正の平均順位和（Rank Sum）


■ t値、z値

\quad t=\dfrac{S}{\sqrt{\dfrac{n \cdot V(T)-S^{2}}{n-1}}}

\quad z=\dfrac{|S| - Correct}{\sqrt{V(T)}}

　※ Correct=連続修正。有：0.5、無：0。



計算例

架空データで計算。

① A列とB列が観測値。C列が前後差。

② 前後差の絶対値に平均順位を付ける。Wilcoxon法はD～I列、Pratt法はJ～O列。
　 ・abs(diff)：前後差の絶対値
　　　D列=IF(C5<>0,ABS(C5),"")
　　　J列=ABS(C5)
　 ・RANK.AVG：平均順位
　　　E列=IF(C5<>"",RANK.AVG(D5,D:D,1),"")
　　　K列=IF(C5<>"",RANK.AVG(J5,J:J,1),"")
　 ・Positive RANK.AVG：前後差が正の平均順位
　　　F列=IF(C5>0,RANK.AVG(D5,D:D,1),"")
　　　L列=IF(C5>0,RANK.AVG(J5,J:J,1),"")
　 ・Negative RANK.AVG：前後差が負の平均順位
　　　G列=IF(C5<0,RANK.AVG(D5,D:D,1),"")
　　　M列=IF(C5<0,RANK.AVG(J5,J:J,1),"")
　 ・t：タイの数
　　　H列=IF(D5<>"",IF(COUNTIF($E$4:E5, E5)=1, COUNTIF(E:E,E5),""),"")
　　　N列=IF(COUNTIF($K$4:K5, K5)=1, COUNTIF(K:K,K5),"")
　 ・(t3-t)：タイの補正項用
　　　I列=IF(H5<>"", H5^3-H5, "")
　　　O列=IF(N5<>"", N5^3-N5, "")

③ 前後差が正の平均順位和 T^{+} を用いて検定統計量 S を算出し、p値を得る。
　 今回は前後差ゼロを３つ含めたので各手法で結果が異なる。


Rで検算。
coinパッケージのデフォルトは、Pratt法・Z近似のp値を返す。連続修正はなし。
zero.method="Wilcoxon" でWilcoxon法。distribution="exact"でexact。

> library(coin)
> ads0 <- data.frame(
+   before=c(10,13,12,12,20,19,18,13,16,20,10,14,20,17,20,18,20,15,13,20,13,18,17,23,16)
+   ,after=c(17,18,16,15,13,13,18,20,10,19,10,20,16,13,13,18,14,13,17,11,11,17,11,12,12)
+ )
> coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="asymptotic")

    Asymptotic Wilcoxon Signed-Rank Test

data:  y by x (pos, neg) 
     stratified by block
Z = 1.4171, p-value = 0.1564
alternative hypothesis: true mu is not equal to 0

> coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Pratt", distribution="asymptotic") #--default

    Asymptotic Wilcoxon-Pratt Signed-Rank Test

data:  y by x (pos, neg) 
     stratified by block
Z = 1.4988, p-value = 0.1339
alternative hypothesis: true mu is not equal to 0


プログラムコード

library(coin)

#-- Z近似
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="asymptotic")
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Pratt", distribution="asymptotic") #--default
  
#-- Exact
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Wilcoxon", distribution="exact")
coin::wilcoxsign_test(before ~ after, data=ads0, zero.method="Pratt", distribution="exact")



proc univariate data=ads;
    var diff;
    output out=out1 probs=probs;
    ods output TestsForLocation=out2;
run;


SASでは、n<20はexact, n>=20ならWilcoxon法・t近似のp値を返す。
（ただし、JMPはPratt法のよう。）

参考

SAS Help Center
JMP Help
https://support.sas.com/resources/papers/proceedings-archive/SUGI94/Sugi-94-172%20Tian.pdf
Solved: Wilcoxon signed rank test : Test statistic S - JMP User Community
Wilcoxon Test • Simply explained - DATAtab
https://www.sciencedirect.com/topics/mathematics/wilcoxon-signed-rank-test

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。