답변 해 주셔서 감사합니다.
다만, 수정 포아송 회귀는 로버스트 분산을 사용하는 것으로 올바른 신뢰 구간을 계산하고 있는 것이지요?
네. 좀 더 자세하게 「올바른」을 말하면, 「표본 사이즈가 큰 경우에 맞다」라고 말할 수 있다고 생각합니다. 된다”라고 말할 수 있다고 생각합니다.통계학의 저곤을 사용하면, “비록 응답 변수의 확률 분포가 무엇이어도, 샌드위치 분산은, 어느 특정의 모델(이번은 Poisson 분포)을 가정한 모델의 추정량의 분산 공분산 행렬에 대한 일치 추정량이 되고 있다”라고 말할 수 있다고 생각합니다.
견고한 분산과 샌드위치 분산은 같은 것을 가리킨다는 이해로 좋을까요?
예. 적어도 이번 문맥에 한정하면, 「로버스트 분산」과 「샌드위치 분산」으로 나는 정확히 같은 것을 가리키고 있습니다.
■ 말씀드린 대로 수정 Poisson 회귀(modified Poisson regression)는 일단 Poisson 회귀를 최대 우도 추정한 후, 로버스트 분산에 의해 표준 오차를 구하고, 그 샌드위치 표준 오차로부터 검정이나 신뢰 구간을 구하고 있습니다.
Poisson 회귀의 최대 우도 추정에서는 Poisson 분포를 가정하므로 이항 분포가 가정되지 않았습니다.
■ 특히 역학에서 수정 Poisson 회귀는 다음과 같은 동기 부여에서 제안된 방법이라고 생각합니다.
어떤 이벤트가 발생할 확률을 p로 했을 경우, 자주 사용되는 모델에 로지스틱 모델이 있습니다.
x=1일 때의 좌변을 log(p1/(1-p1)), x=0일 때의 좌변을 log(p0/(1-p0))로 하면, β1은,
β1 = log(p1/(1-p1)) - log(p0/(1-p0))
를 나타냅니다.
β1 = log({p1/(1-p1)}/{p0/(1-p0)})
가 됩니다. 양변을 지수 변환하면, exp(log(a))=a 보다,
exp(β1) = {p1/(1-p1)}/{p0/(1-p0)}
즉, β1은 (x가 1 변화했을 때의) 대수 확률비를, exp(β1)는 (x가 1 변화했을 때의) 확률비를 나타내고 있습니다.
{p1/(1-p1)}/{p0/(1-p0)}≒p1/p0
확률 비율은 위험 비율 (비율 비율)에 근사 할 수 있습니다. 그러나이 근사는 p1, p0가 0에 가깝지 않으면 성립하지 않습니다. 르의 우변이 어떤 바람이 되어도 p는 (0 이상 1 이하라고 함) 확률의 성질을 거의 거의 유지합니다 (단, 로지스틱 회귀 모델에서는 확률이 정확히 0이 될 때와 1이 될 때 문제가 발생합니다).
확률비의 해석이 어렵기 때문에 분석자가 해석하기 쉽도록 리스크비에 대응한 모델은
와 같이, 좌변을 로짓트가 아니고, 비율의 대수로 한 것입니다.
(로짓트(log(p/(1-p)))나 로그(log(p)) 등의 변환식을 링크 함수라고 합니다.)
이 때, log(p1) - log(p0) = β1 이 되므로,
log(p1/p0) = β1
p1/p0 = exp(β1)
즉, 이 로그 링크 함수의 모델에서는 β1은 대수 리스크비, exp(β1)가 리스크비가 됩니다.
그러나 이 모델은 문제가 있습니다.우변이 예를 들어 β0 + β1x = 1.0이 되면, log(p)=1.0 즉 p = exp(1.0)≒2.72가 되어, 확률인데 p가 1을 넘는 때가 있습니다.
이 문제를 피하기 위해 링크 함수가 로그 인 경우와 궁합이 좋은 Poisson 분포를 기반으로 최대 우도 추정을 할 수 있습니다. 그러나 그렇게 하면 (본래는 이항 분포를 따르고 있다고 생각되는데) Poisson 분포에서 표준 오차·검정·신뢰 구간이 계산되게 됩니다. 따라서 응답 변수의 확률 분포 지정에 견고한 샌드위치 분산 (견고한 분산 및 경험 분산)을 사용하는 것이 제안되었습니다. 샌드위치 분산을 지정하면 응답 변수의 확률 분포로 잘못된 것을 지정해도 매개변수 추정량에 대한 샌드위치 표준 오차(및 샌드위치 공분산 행렬)는 표본 크기가 크면 합리적이 됩니다(참 표준 오차에 대한 일치 추정량이 됩니다).
덧붙여 「확률비의 해석이 어렵다」라고 하는 것은 아마 2점 있다고 생각합니다. (1-b))」라고 직감적으로는 무엇이 무엇인지 조금 모릅니다.또, 리스크 차이나 리스크비는 병합 가능성(collapsibility)이라고 하는 성질이 있는 것에 대해, 오즈비는 병합 가능성이 성립하지 않습니다.
상기와 같이, 수정 Poisson 회귀가 주목받고 있는 것은, 「승률비라고 해석이 어렵다」라고 하는 단점을 보충하는 것이므로, 확률비의 해석이 별로 어렵지 않으면, 종래의 로지스틱 회귀 모델인 채로 좋다 는 없을까 개인적으로는 생각합니다.(예를 들어, 분석의 목적이 「예측」이나 「분류」이고, 회귀 계수의 해석이 중시되지 않는 경우, 로지스틱 회귀 쪽이 수정 Poisson 회귀보다 사용하기 편하다고 생각합니다.)
■전날, 말한 방법(R이나 Python을 호출하는, 및,JMP스크립트 언어로 스크래치로 프로그램 한다)는, 프로그래밍이 필요하므로, JMP의 가장 좋은 점(노코드로 대화적으로 즐겁게 통계 분석을 할 수 있다)가 사라져 버립니다.그러나, JMP의 유저라도, 특히 루틴적인 통계 처리를 실시하고 있는 분으로.
JMP Pro 19가 되면 마우스 조작으로 수정 Poisson 회귀를 할 수 있게 됩니다.현재로 JMP상에서 실시하고 싶은 경우에는 아무래도 프로그래밍이 필요합니다.
이하에, 작은 의사 데이터를 샘플 데이터로서,
① R을 JMP에서 호출하여 수정 Poisson 회귀를 실행하는 방법
② Python을 JMP에서 호출하여 수정 Poisson 회귀를 실행하는 방법
③JMP 스크립트 언어로 스스로 제로로부터 프로그래밍해 수정 Poisson 회귀를 실행하는 방법
우선은, 대략 봐 주셔, 어느 정도 귀찮은 것의 대략적인 기준을 체감해 주셨으면 합니다.아래의 코드는, 어디까지나 개인이 작성한 하나의 샘플이며, 폐사가 조직으로서 제공하는 제품은 전혀 없습니다(책임은 무엇입니까)
게다가 아래의 설명은 여기에서 다루는 작은 의사 데이터만으로 움직이는 설명입니다.
여기에서는 설명을 위해 상당히 작은 데이터를 사용하고 있습니다만, 본래, 샌드위치 분산은 어느 정도의 대규모 데이터에 대해서 합리적인 방법이 되고 있습니다(점근론에 근거하는 방법이 되고 있습니다).
①, ②에 대해서는, 열명에 일본어가 들어가 있으면 잘 움직이지 않을지도 모릅니다.또, 나 자신, R나 Python에 자세하지 않기 때문에, 보다 우아한 패키지나 프로그래밍 방법이 있을지도 모릅니다.
③에 대해서는 다음 사항을 전혀 고려하지 않습니다.
a) 결측값 처리: 데이터에 결측값이 있으면 오류가 발생합니다.
b) 계획 행렬의 랭크 낙하에 대한 에러 처리: 선형 회귀 분석이나 Poisson 회귀를 비롯한 일반화 선형 모델에서는 설명 변수 사이에 1차 종속성이 있으면 매개변수 추정치가 고유하게 정해지지 않습니다.
c) 분리에 대한 에러 처리: 로지스틱 회귀나 Poisson 회귀에는 데이터에 따라서는 분리(separation)라는 문제가 생겨 파라미터 추정값이 +∞ 또는 -∞가 되는 경우가 있습니다.
①과 ③을 하나로 정리한 JSL 코드의 파일을 첨부해 둡니다.
②의 Python 코드의 파일도 첨부해 둡니다.이것을 실행하려면 , 상기의 JSL 파일의 1~7행 부분의 코드로 JMP 데이터 테이블이 열려 있지 않으면 안됩니다.
■0) 우선, 이하의 ①, ②, ③에서 이용하는 샘플 데이터를 작성합니다.
0-1) 파일(F) > 새로 만들기(N) > JSL 스크립트(Cntl + T)를 선택합니다.
0-2) 호출 된 스크립트 편집기에서 마우스 오른쪽 버튼을 클릭하고 팝업 메뉴에서 [창에 로그 표시]를 선택합니다.
0-3) 스크립트 편집기의 상단에 다음 코드를 복사합니다.이 코드는 샘플 데이터를 만듭니다.
New Table("Test",
New Column("X1", Nominal, Character,Values({"A","A","A","A","A","A","B","B","B","B","B","B","B","B","B"})),
New Column("X2", Values([5,4,3,2,1,4,5,4,3,3,3,2,1,5,4])),
New Column("Y", Values([0,1,0,1,0,1,1,0,1,1,0,1,1,0,1]));
);
0-4) 스크립트 편집기에서 마우스 오른쪽 버튼을 클릭하고 팝업 메뉴에서 스크립트 실행을 선택합니다.
0-5) 위의 데이터 값을 가진 JMP 데이터 테이블이 생성됩니다.
■①R을 JMP로부터 호출해 수정 Poisson 회귀를 실행하는 방법
1) R 설치 및 JMP에서 R이 호출되는지 확인
먼저 JMP에서 R을 호출하려면 JMP (JMP, JMP Pro, JMP Student Edition)가 설치된 시스템에 R이 설치되어 있어야합니다.
R이 벌써 인스톨 되고 있는 경우, 혹은, R를 인스톨 한 후에, JMP로부터 R가 불려 갈지 어떨지를 체크합니다.
1-1) 다음 코드를 스크립트 편집기에 복사합니다.
R Init();
R Submit("
R.Version();
print('Hello, World!')
");
1-2) 스크립트 편집기에서 마우스 오른쪽 버튼을 클릭하고 팝업 메뉴에서 스크립트 실행을 선택합니다.
1-3) 스크립트 에디터 아래의 로그에 R의 버전이나 "Hello, World!"라는 문자열이 나오면 JMP에서 R을 호출할 수 있습니다.
2) 다음으로 JMP 데이터 테이블을 위의 R Init() 에 의해 JMP 와 연결된 R 에 보냅니다.
2-1) 데이터 테이블을 연 상태에서 다음 코드를 실행합니다 (이미 R에 연결되어 있으면 첫 번째 R Init() 필요 없음).
R Init();
dt = Current Data Table();
dt << New Column("ID", Nominal, Formula(Row()));
Wait(0);
R Send(dt);
Current Data Table()은 현재 현재 데이터가 되고 있는 JMP 데이터 테이블에 대한 참조를 가져옵니다.
dt << New Column("ID", Nominal, Formula(Row()));"은 각 행에 대한 일련 번호의 열을 만듭니다. 이미 일련 번호(ID 번호)의 열이 있는 경우에는 이 한 줄이 필요하지 않습니다.
마지막 R Send(dt)에서 데이터 테이블을 R로 보냅니다.
3) R 코드를 실행합니다.
R에 데이터를 보낼 수 있으면 나중에 R에서 수정 Poisson 회귀를 추정합니다.
GEE(일반화 방정식 추정)는 원래 도수 데이터 등에서 상관이 있는 경우의 추정 방법으로서 제안된 것입니다만, 거기서 표준 오차를 계산할 때에 샌드위치 분산이 사용되고 있었습니다. 하는 기능을 이용할 수 있습니다.이것을 위해서는, ID열로서는 각 행에 일의인 것을 이용합니다 (또, 상관 구조에는 독립을 이용합니다).
아래의 코드에서는 geeglm을 이용하여 X1과 X2를 설명변수로 한 수정 Poisson 회귀를 하고 있습니다(GEE의 말로 말하면 로그 링크 함수로 Poisson 분포를 가정한 모델로 각 행을 독립으로 간주하여 샌드위치 분산을 구했다).
결과는 스크립트 창 하단의 로그에 출력됩니다.
R Submit("
if (!require('geepack')) {
install.packages('geepack')
}
library(geepack)
gee.result <- geeglm(Y~X1 + X2, data=dt, id = ID, family=poisson('log'))
summary(gee.result)
est <- coef(gee.result)
se <- diag(vcov(gee.result))
z <- qnorm(0.975)
print('Risk Ratios and the Confidence Interlvals')
print(cbind(exp(est), exp(est - se*z), exp(est + se*z)))
");
다음은 rqlm 패키지를 사용한 예입니다. 이것은 다음의 외부 자료를 참고하여 프로그래밍한 것입니다.
노마 히사시(2025)
rqlm : 수정 포아송 회귀 및 수정 최소 제곱 회귀 분석을위한 R 패키지
/* 以下のURLを参考にしました。https://www.ism.ac.jp/~noma/file/software/rqlm_JPN.pdf */
R Submit("
if (!require('rqlm')) {
install.packages('rqlm')
}
library(rqlm)
mpr.result <- rqlm(Y ~ X1 + X2, data=dt, family=poisson, eform=TRUE)
mpr.result
");
■①Python을 JMP로부터 호출해 수정 Poisson 회귀를 실행하는 방법
JMP18부터 JMP를 설치하면 JMP 내부에서만 작동하는 Python이 함께 설치됩니다. 따라서 Python을 새로 설치할 필요가 없습니다.
1-0) 이하에서는 Python의 코드로 프로그래밍해 나가지만, JMP상에서 Python 코드를 작성하려면 , 파일(F) > 신규 작성(N) > Python 코드를 선택합니다.
1-1) 몇 가지 필요한 Python 패키지를 설치합니다. 이번에 필요한 패키지는 numpy, pandas, statsmodels입니다.
from jmputils import jpip
jpip('install --upgrade', 'pip setuptools certifi')
jpip('install', 'numpy pandas statsmodels')
이 코드는 한 번 실행하면 다음 번부터 실행할 필요가 없습니다.
1-2) 패키지와 함수를 가져옵니다.
다음 코드를 Python 스크립트 편집기에 복사하여 실행하십시오.
import jmp
import numpy as np
import pandas as pd
import statsmodels as sm
from statsmodels.genmod.generalized_estimating_equations import GEE
from statsmodels.genmod.families import Poisson
from statsmodels.genmod.families.links import log
from statsmodels.genmod.cov_struct import Independence
1-3)
JMP 데이터 테이블을 Python으로 읽어 들여서 pandas의 데이터 프레임으로 변환합니다.
이 예에서는 현재 (현재) JMP 데이터 테이블 (지금 열려있는 JMP 데이터 테이블)을 읽고 데이터를 pandas의 데이터 프레임으로 만듭니다.
데이터에 따라 다음 코드를 변경해야 합니다.
dt = jmp.current()
data = pd.DataFrame(
{
'Y': dt['Y'][0:n],
'X1': dt['X1'][0:n],
'X2': dt['X2'][0:n],
'ID': dt['ID'][0:n]
}
)
파이썬 GEE 함수를 사용하려면 직접 지시 변수 (더미 변수)를 만들어야합니다.
이 예에서는 X1이 이진 변수입니다. 여기서 X1 == "B"일 때 1, 그렇지 않으면 0 인 더미 변수를 만듭니다.
data['X1B'] = (data['X1'] == 'B').astype(int)
print(data)
1-4)
GEE를 실행합니다. 응답 변수의 확률 분포에 Poisson 분포를, 링크 함수에 로그를 지정합니다.
또한 파라미터 추정값과 그 신뢰 구간을 지수 변환하고 있습니다.
model = GEE.from_formula( "Y ~ X1B + X2", groups="ID",
data=data,
family=Poisson(link=log()),
cov_struct=Independence())
model = GEE(Ydata, Xdata, groups=ID, cov_struct=cov_struct, family=Poisson(link=Log()))
result = model.fit()
print(result.summary())
print(result.cov_robust)
print(np.exp(result.params))
print(np.exp(result.conf_int()))
■③ JMP 스크립트 언어로 스스로 제로로부터 프로그래밍해 수정 Poisson 회귀를 실행하는 방법
수정 Poisson 회귀는 먼저 Poisson 회귀의 최대 우도 추정에서 점 추정치를 찾습니다.
이때는 Newton-Raphson법이라는 방법이 자주 사용되고 있습니다. 모델))이라는 모델족에 속하고 있지만 일반화 선형 모델에서는 이 기울기와 헤세 행렬의 계산이 간단해집니다.
또, 스스로 제로로부터 프로그램 하는 경우에는, 명의 척도의 변수는 스스로 지시 변수(더미 변수)를 작성할 필요가 있습니다.
세세한 점의 설명은 생략합니다만, 에러 처리(결측치를 어떻게 하는지, 랭크 낙하를 어떻게 하는지, 분리를 어떻게 하는지)를 무시하면, 이하 정도의 코드량으로 수정 Poisson 회귀는 실시할 수 있습니다.
Names Default To Here(1);
dt = Current Data Table();
/* 自分で指示変数（ダミー変数）を作成しないといけない */
x1 = Design(Column(dt, "X1")<
<GET as="" matrix=""></GET>
<DIV>설명이 길어져서 죄송합니다. Marketplace라고 하는 사이트에서 제공되고 있습니다(여기에서 제공되고 있는 확장 기능은, 현재 상태 모습(as-is)으로 제공되고 있는 것입니다) 유감스럽지만, 수정 Poisson 회귀의 기능은 JMP Marketplace에서도 제공되고 있지 않습니다.</DIV>
<DIV>
</DIV>
<DIV>JMP인 처리가 가능한지, 할 수 없는가에 대해서, 이하와 같은 레벨로 나눌 수가 있을지도 모릅니다. 이번은, 레벨 d)에 상당하는 것이라고 생각합니다.</DIV>
<DIV>a) JMP나 JMP Pro로 보통 준비되어 있는 기능 → 순식간에 할 수 있다</DIV>
<DIV>
<DIV>b) JMP나 JMP Pro의 기능을 이용하는 경우 → 궁리하면 바로 할 수 있다</DIV>
</DIV>
<DIV>
<DIV>c) JMP Marketplace에서 제공되는 기능 → 어디까지나 현재 상태 모습이지만 바로 할 수 있다</DIV>
<DIV>d) 위에서 제공되지 않은 기능으로 약간 간단한 것 → 프로그래밍하면 할 수 있다 (특히 R이나 Python으로 할 수 있으면)</DIV>
<DIV>e) 상기에서 제공되지 않은 기능으로 상당히 고도의 것 → 상당한 전문 지식이 있는 경우에만 가능</DIV>
<DIV>
</DIV>
<DIV>마지막으로 일점만, 신경이 쓰였으므로, 서두의 「정확함」에 대해서 보충해 둡니다.원래, 회귀 분석 자체, 제일 간단한 최소 제곱의 선형 회귀도 포함해, 여러가지 사용되는 방법이 있습니다.</DIV>
<DIV>이와사키학(2021)</DIV>
<DIV>통계적 인과 추론의 관점에 따른 중회귀 분석</DIV>
<DIV>일본 통계학회지, 50(2), 363-379</DIV>
<DIV>
<A class="bluelink-style" href="https://doi.org/10.11329/jjssj.50.363" target="_blank" rel="nofollow noopener noreferrer">https://doi.org/10.11329/jjssj.50.363</A>
</DIV>
<DIV>위의 여기까지의 설명은 어디까지나 수정 Poisson 회귀의 실행 방법을 말했을 뿐, 지금 놓인 상황에서, 회귀 분석을 실시하는 것이 「올바른」인가, 수정 Poisson 회귀를 실시하는 것 가 '올바른'인지, 게다가 현재 상황이 회귀분석이 가정하고 있는 큰 전제를 충족하고 있다는 의미에서 '올바른'인지 등의 의미에서 '정확성'을 실현하는 것은 전혀 없습니다. 현실의 상황에 대해 현재의 통계 분석이 「올바른」인가는, (나 자신은 판단하는 것이 전혀 할 수 없기 때문에) 그 분야의 응용 통계가에게 상담하는 것을 추천하고 싶습니다 (나 자신, 컨설팅을 받은 적은 없지만, 분석 그 자체의 대안을 제시하는 일도 있을까요, 또, 수정 Poisson 회</DIV>
<DIV>
</DIV>
</DIV>
