17/12/21 05:20:38.89 NCbCbV7K.net
> 0.05^(1/(1:100+1))
[1] 0.2236068 0.3684031 0.4728708 0.5492803 0.6069622 0.6518363 0.6876560
[8] 0.7168712 0.7411344 0.7615958 0.7790778 0.7941833 0.8073638 0.8189637
[15] 0.8292503 0.8384339 0.8466824 0.8541315 0.8608917 0.8670541 0.8726946
[22] 0.8778766 0.8826538 0.8870719 0.8911696 0.8949808 0.8985343 0.9018554
[29] 0.9049661 0.9078859 0.9106318 0.9132188 0.9156604 0.9179684 0.9201535
[36] 0.9222253 0.9241923 0.9260624 0.9278425 0.9295389 0.9311574 0.9327032
[43] 0.9341812 0.9355957 0.9369507 0.9382499 0.9394966 0.9406940 0.9418449
[50] 0.9429520 0.9440178 0.9450445 0.9460342 0.9469889 0.9479105 0.9488005
[57] 0.9496607 0.9504924 0.9512971 0.9520760 0.9528305 0.9535615 0.9542703
[64] 0.9549577 0.9556248 0.9562724 0.9569014 0.9575126 0.9581067 0.9586843
[71] 0.9592463 0.9597932 0.9603256 0.9608441 0.9613492 0.9618415 0.9623214
[78] 0.9627893 0.9632458 0.9636912 0.9641260 0.9645504 0.9649650 0.9653699
[85] 0.9657656 0.9661524 0.9665305 0.9669003 0.9672620 0.9676158 0.9679621
[92] 0.9683011 0.9686330 0.9689580 0.9692763 0.9695882 0.9698938 0.9701933
[99] 0.9704870 0.9707748
>
337:卵の名無しさん
17/12/22 02:11:19.86 J9UAx7pH.net
シンプソンのパラドックス
ある仮想疾患の治癒率
軽症 重症
国立大学 10/10 10/90
底辺私立 70/90 0/10
自然経過 40/50 5/50
国立大学の方が軽症・重症とも成績がよいが
総数比較では底辺私立の方が成績がよい。
この疾患は自然治癒率が45%とされています。
この疾患の底辺私立での治癒率は70%です。
これに�
338:ﾎして国立大学での治癒率はわずか20%です。 という記述も嘘ではないね
339:卵の名無しさん
17/12/22 12:38:20.39 FPEZRkaT.net
f <- function(n=10,alpha=1,beta=1,Print=FALSE){
N=n
z=n
if(Print) {
bayes=binom::binom.bayes(z,N, prior.shape1=alpha,prior.shape2=beta)
show(binom::binom.bayes.densityplot(bayes))
}
hdi=HDInterval::hdi(qbeta, shape1=alpha+z, shape2=beta+N-z)
return (c(lower=hdi[1],mean=(alpha+z)/(alpha+beta+N),
mode=(alpha+z-1)/(alpha+beta+N-2), upper=hdi[2]))
}
f(10,P=TRUE)
nn=1:30
yy=sapply(nn,function(x)f(x,Print=FALSE)[1])
plot(nn,yy,pch=19,xlab='裏口バカ連続合格数',ylab='裏口確率信頼区間下限')
curve((0.05)^(1/(x+1)),add=TRUE,lty=3) # 0.05の（合格者数＋１）乗根
340:卵の名無しさん
17/12/22 20:41:25.07 JBU22EfC.net
# N回続けて裏、事前分布はmode値0.5, 集中度(形状母数和)=kappa
source('tools.R')
N=5 ; z=5
Bayes(N,z,alpha=1,beta=1,Print=TRUE) ; 0.05^(1/(N+1)) # N=z
# 事前分布が最頻値0.5で集中度(κ=α+β)のとき事後分布の関係
.mode=0.5
Kappa2Bayes <- function(kappa,.mode=0.5){
AB=betaABfromModeKappa(.mode,kappa)
Bayes(N,z,alpha=AB[[1]],beta=AB[[2]])
}
K=seq(2,1000,by=0.5)
K=K[-1]
res=sapply(K,Kappa2Bayes)
Mat=as.matrix(res)
plot(K,Mat['lower',],type='l',xlab=bquote(kappa),
ylab='Probability', ylim=c(0,1),lty=3)
lines(K,Mat['mean',],col=4,lwd=4)
lines(K,Mat['mode',],col=2,lwd=2)
lines(K,Mat['upper',],lty=3)
legend('bottomright',bty='n',legend=c('mean','mode','upr','lwr'),
col=c(4,2,1,1),lty=c(1,1,3,3),lwd=c(4,2,1,1))
341:卵の名無しさん
17/12/22 21:31:09.94 JBU22EfC.net
N(=100)回コインをなげてｎ(=5回)以上続けて表がでる
seqn<-function(n=5,N=100,p=0.5){
rn=rbinom(N,1,p)
count=0
for(i in 1:N){
if(rn[i] & count<n){
count=count+1
}
else{
if(count==n) {return(TRUE)}
else{
count=0
}
}
}
return(count==n)
}
mean(replicate(10^6,seqn()))
> mean(replicate(10^6,seqn()))
[1] 0.810223
342:卵の名無しさん
17/12/23 05:59:24.55 RXgTWS3Z.net
pooledVariance <- function(...) {
args = list(...)
n.args=length(args)
ss2=0
df=0
for(i in 1:n.args){
ss2 = ss2 + var(args[[i]])*(length(args[[i]])-1)
df = df + (length(args[[i]])-1)
}
ss2/df
}
effectsize <- function(y1,y2){
diff=mean(y1)-mean(y2)
var=(var(x1)*(length(x1)-1)+ var(x2)*(length(x2)-1))/(length(c(y1,y2))-2)
sd=sqrt(var)
diff/sd
}
library(effsize)
cohen.d()
343:卵の名無しさん
17/12/23 06:36:30.75 RXgTWS3Z.net
solve[{M=(a-1)/(a+b-2), V=a*b/((a+b)^2*(a+b+1))},{a,b}]
344:卵の名無しさん
17/12/23 17:39:59.62 K1KZza+Z.net
> 1/(1-(1-0.99)^(1/317))
[1] 69.33689
1 / ( 1- n√(1-confidence.level) )
345:卵の名無しさん
17/12/23 18:01:15.54 K1KZza+Z.net
# 1 / ( 1- n√(1-confidence.level) )
confidence.level=0.95
rule3 <- function(n,confidence.level=0.95){ # n人に1人の副作用
p=1/n
q=1-p # q^n.sample < 1-confidence.level
n.sample = log(1-confidence.level)/log(q)
return(n.sample)
}
nn=seq(1,10000,by=100)
plot(nn,sapply(nn,rule3))
curve(3*x,col=2,add=TRUE)
plot(nn,sapply(nn,function(x) rule3(x,conf=0.99)))
lm( sapply(nn,function(x) rule3(x,conf=0.99)) ~ nn + 0)
curve(4.605*x,col=4,add=TRUE)
346:卵の名無しさん
17/12/24 20:11:33.25 SYFul+nD.net
サイコロの１の目が何回続けてでたらイカサマサイコロかを考えて暇つぶし。
このスレの趣旨ｗに合わせてこんな問題にしてみた。
ド底辺学力学生がド底辺特殊シリツ医大を受験したとする。
試験は五者択一で三教科150問、合格ラインは６割とされる。
中学数学すらできないド底辺学力ゆえ正解できるのは偶然に頼るしか�
347:ﾈく、 正答率は概ね1/５で、その日の運で1/4から1/6と推定されている。 これを、正答確率は最頻値1/５で1/6から1/4の間に正答する確率の95%があると設定する。 このド底辺が合格したとする。150問中90問以上正答したことになる。 これがイカサマ入試である確率を求めよ。 事前確率のβ分布のパラメータ算出がやや手間だが、あとはルーティン作業 JAGSを使ってMCMCでの結果 http://i.imgur.com/V7TaBG7.png 解析解： 0.9994608 とほぼ一致。
348:卵の名無しさん
17/12/25 08:52:06.92 Nj//P1mP.net
require(rjags)
N=10
z=0
y=c(rep(1,z),rep(0,N-z))
ph=c(1,2/3,1/2,1/5)
pc=c(2/100,50/100,40/100,8/100)
names(ph)=c('特待生','学力合格','加点合格','ガチの裏')
names(pc)=c('特待生','学力合格','加点合格','ガチの裏')
dataList5=list(N=N,y=y,ph=ph,pc=pc)
# JAGS model
modelString5 ="
model {
for(i in 1:N){
y[i] ~ dbern(ph[coin])
}
coin ~ dcat(pc[])
}
"
writeLines(modelString5,'TEMPmodel.txt')
jagsModel5=jags.model('TEMPmodel.txt',data=dataList5)
codaSamples5=coda.samples(jagsModel5,var=c('coin'),n.iter=100000,na.rm=TRUE)
summary(codaSamples5)
js5=as.matrix(codaSamples5)
re=numeric()
for(i in 1:4) re[i]=mean(js5==i)
dat=data.frame(割合=round(re*100,3))
rownames(dat)=names(pc)
dat
349:卵の名無しさん
17/12/25 16:31:55.67 UMwuImpO.net
頻度主義統計の謎。
立方体からなるサイコロの目のでる確率はすべて等しく1/6である、を帰無仮説とする。
そのサイコロをふって1の目がでた。2回目は2の目がでた。
その確率は1/6*1/6で1/36=0.02778 < 0.05だから帰無仮説は棄却される。
どの目の組合せでも同じく帰無仮説は棄却される。
頻度主義統計のもとではすべてのサイコロはいびつである。
350:卵の名無しさん
17/12/27 10:24:07.92 un/eaZi1.net
a005=0.05
n100=100
HDInterval::hdi(qbeta,shape1=n100+.a,shape2=.b)
qbeta(a005,n100+.a,.b)
a005=0.05
# p^n100 < a005
# p < a005^(1/n100)
GolgoLowerLimit <- function(a005,n100=100){ # Golgo lower limit
c(a005^(1/(n100+1)),qbeta(a005,n100+1,1))
}
GolgoLowerLimit(0.05)
AHO <- function(a00,n100=100,shoot=10000){
(a005^(1/n100)+1)/2*shoot
}
AHO(0.05)
x=seq(0.001,0.1,by=0.001)
plot(x,sapply(x,function(x) AHO(x,100,10000)),type='l',lwd=2,
las=1,ylab='「平均値」',xlab='危険率')
351:卵の名無しさん
17/12/27 10:57:46.05 un/eaZi1.net
data{
int Npip; // 6
real alpha[Npip]; // c(1,1,1,1,1,1)
int Ntotal; // length(y)
int y[Ntotal];
}
parameters{
simplex[Npip] pi;
}
model{
for(i in 1:Ntotal){
y[i] ~ categorical(pi);
pi ~ dirichlet(alpha)
}
}
352:卵の名無しさん
17/12/28 15:27:03.23 t4TEzXKz.net
source('tools.R')
RBI <- function(a=1,b=1,zi=1,Ni=1,zj=0,Nj=2,ROPE=NULL,Print=TRUE){
(y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj)))
(s=as.numeric(factor(c(rep('D',Ni),rep('U',Nj)))))
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
# JAGS model
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( sIdx in 1:Nsubj ) {
theta[sIdx] ~ dbeta(", a,',' , b," )
}
}
")
# close modelString
353:卵の名無しさん
17/12/28 15:41:40.10 t4TEzXKz.net
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
mcmcMat=as.matrix(
354:codaSamples) if(Print){ par(mfrow=c(2,1)) curve(dbeta(x,a,b),xlab=paste0('Beta(',a,',',b,')'),bty='n',yaxt='n',ylab='', type='h', col='blue') print(plotPost2(mcmcMat[,1]-mcmcMat[,2],compVal=0,ROPE=ROPE, cenT='mean',xlab=bquote(Delta),cex=1,showCurve=FALSE)) } dif=mcmcMat[,1]-mcmcMat[,2] print(c(HDInterval::hdi(dif),mean=mean(dif))) invisible(dif) }
355:卵の名無しさん
17/12/29 18:23:01.86 lWMSqtax.net
URLﾘﾝｸ(image.slidesharecdn.com)
356:卵の名無しさん
17/12/29 18:23:24.53 lWMSqtax.net
URLﾘﾝｸ(image.slidesharecdn.com)
357:卵の名無しさん
17/12/30 02:13:47.38 8e/jZDFc.net
BF01<- function(n,z,p0,a=1,b=1,p1=0.5) {
bf=gamma(a)*gamma(b)/gamma(a+b) *
gamma(a+b+n)/gamma(a+z)/gamma(b+n-z) *
p0^z*(1-p0)^(n-z)
pri=p1/(1-p1)
pos=pri*bf
post.prob=pos/(1+pos)
c(BayesFactor=bf, PostProb=post.prob)
}
BF01(5,4,0.5,1,1)
BF01(7,7,0.5,1,1)
358:卵の名無しさん
17/12/31 08:13:46.23 2OFZ1/Lf.net
URLﾘﾝｸ(to-kei.net)
359:卵の名無しさん
17/12/31 08:29:05.13 2OFZ1/Lf.net
a=1
b=1
n=5
z=5
p=0.5
dataList=list(
a=a,
b=b,
n=n,
z=z,
p=p
)
modelString,='
model{
bf=gamma(a)*gamma(b)/gamma(a+b) *
gamma(a+b+n)/gamma(a+z)/gamma(b+n-z) *
p^z*(1-p)^(n-z) # bf= (n+1)*p^n
pri=p0/(1-p0)
pos=pri*bf
post_prob=pos/(1+pos)
p0 ~ beta(1,1)
}
'
360:卵の名無しさん
17/12/31 09:49:21.42 2OFZ1/Lf.net
a=1
b=1
n=5
z=n
p=0.5
dataList=list(
a=a,
b=b,
n=n,
z=z,
p=p
)
modelString,='
model{
bf= (n+1)*p^n
pri=p0/(1-p0)
pos=pri*bf
post_prob=pos/(1+pos)
p0 ~ beta(1,1)
}
'
writeLines(modelString,'TEMPmodel.txt')
jagsModel=jags.model('TEMPmodel.txt', data=dataList)
codaSamples=coda.samples(jagsModel,n.iter=10000,var=c
('p0'))
js=as.matrix(codaSamples)
361:卵の名無しさん
17/12/31 10:18:36.31 2OFZ1/Lf.net
５試合連続で勝敗予想的中なら頻度主義では予知能力あるとされる。p=0.03125 < 0.05
URLﾘﾝｸ(to-kei.net)
ベイズでやってみるなら
的中率が1/2である確率は一応分布に従う（事前分布）として
５試合連続的中した後の的中率事後分布がどうなるかを考える。
n=5
k=10^4
p0=runif(k)
bf= (n+1)*p^n
pri=p0/(1-p0)
pos=pri*bf
post_prob=pos/(1+pos)
quantile(post_prob,prob=c(0.025,0.5,0.975))
362:卵の名無しさん
17/12/31 16:32:20.16 dbAzKAtn.net
dataList=list(
a=a,
b=b,
n=n,
z=z,
p=p
)
modelString='
model{
bf = p^z*(1-p)^(n-z)*exp(loggam(a))*exp(loggam(b))/exp(loggam(a+b))*exp(loggam(a+b+n))/exp(loggam(a+z))/exp(loggam(b+n-z))
pri = p0/(1-p0)
pos = pri*bf
post_prob = pos/(1+pos)
p0 ~ dbeta(a,b)
}
'
writeLines(modelString,'TEMPmodel.txt')
jagsModel=jags.model('TEMPmodel.txt', data=dataList)
codaSamples=coda.samples(jagsModel,n.iter=20000,chains=4,var=c('post_prob'))
js=as.matrix(codaSamples)
xlim=c(0,min(1,mean(js[,'post_prob'])*6))
BEST::plotPost(js[,'post_prob'],xlim=xlim,xlab=paste(n,'guess',z,'correct'),cenTend='mean',
showCurve = FALSE,col='gray')
363:卵の名無しさん
17/12/31 17:30:53.02 dbAzKAtn.net
# X[n] ~ N(μ,σ)
# prio.μ ~ N(η,τ)
# post.μ ~ N( [(n/σ2)x_ + (1/τ2)η] /[(n/σ2) + (1/τ2)],
# (σ2/n)*τ2/[(σ2/n) + τ2]
post_norm <- function(eta,tau,n,x_,sigma){
mean=(n/sigma^2*x_ + 1/tau^2*eta)/(n/sigma^2 + 1/tau^2)
sd=(sigma^2/n)*tau^2 / ((sigma^2/n) + tau^2)
return(c(mean,sd))
}
post_norm(eta=180,tau=15,n=5,x_=195,sigma=10)
364:卵の名無しさん
17/12/31 17:31:27.73 dbAzKAtn.net
n=5
X_=195
sigma=10
x=rnorm(n,X_,sigma)
X=(x-mean(x))/sd(x)*sigma + X_
mean(X);sd(X)
eta=180
tau=15
dataList=list(n=n,sigma=sigma,X=X,eta=eta,tau=tau)
modelString='
model{
for(i in 1:n){
X[i] ~ dnorm(mu,1/sigma^2)
}
mu ~ dnorm(eta,1/tau^2)
}
'
writeLines(modelString,'TEMPmodel.txt')
jagsModel=jags.model('TEMPmodel.txt',data=dataList)
codaSamples=coda.samples(jagsModel,n.iter = 100000,var=c('mu'))
js=as.matrix(codaSamples)
hist(js)
mean(js)
var(js)
post_norm(eta=180,tau=15,n=5,x_=195,sigma=10)
365:卵の名無しさん
18/01/01 21:30:04.10 9r6fKQha.net
URLﾘﾝｸ(i.imgur.com)
366:卵の名無しさん
18/01/03 20:47:49.71 CFofwzsi.net
■3囚人問題（英: Thre
367:e Prisoners problem） ある監獄にA、B、Cという3人の囚人がいます 3人のうちランダムに選ばれた1人に恩赦が出ます 誰が恩赦になるかは看守は答えない 囚人Aに看守が「Bは死刑になる」と教えてくれます この時、看守は嘘は言いません 囚人Aに恩赦が与えられる確率は何％でしょうか？
368:卵の名無しさん
18/01/03 20:49:01.81 CFofwzsi.net
死刑囚A,B,CでAが看守に尋ねてBは死刑執行されると告げられたと設定。
恩赦(onsha)を受けるをo,死刑執行されると告(tsuge)げられるをtで表す。
Aが恩赦を受ける確率P(A=o)=1/3
Bが恩赦を受ける確率P(B=o)=1/3
Cが恩赦を受ける確率P(C=o)=1/3
求めたいのは、Bが死刑執行されると告げられた後のAが恩赦を受ける確率P(A=o|B=t)である。
ベイズの公式により
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o)
P(B=t|B=o)=0 Bが恩赦を受けるときBが死刑執行されると告げられる確率=0
P(B=t|C=o)=1 CBが恩赦を受けるときBが死刑執行されると告げられる確率=1
問題は P(B=t|A=o)
恩赦を受けるのがAであるときに看守がCではなくBが死刑執行されると告げる確率は示されていない。
この確率をｐとすると
P(A=o|B=t)は p/(p+1)となる。
もちろんp=0.5であれば、P(A=o|B=t)=1/3と看守に告げられる前と同じである。
ここでｐが一様分布からさまざなβ分布に従うとするとどうなるか、グラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
左の緑が看守がBとCが死刑執行予定であるときにBを選んで答える確率分布。
右の青が看守がBと告げたときのAが恩赦を受ける確率の分布。
369:卵の名無しさん
18/01/03 21:28:14.93 CFofwzsi.net
（タイプミス修正）
死刑囚A,B,CでAが看守に尋ねてBは死刑執行されると告げられたと設定。
恩赦(onsha)を受けるをo,死刑執行されると告(tsuge)げられるをtで表す。
Aが恩赦を受ける確率P(A=o)=1/3
Bが恩赦を受ける確率P(B=o)=1/3
Cが恩赦を受ける確率P(C=o)=1/3
求めたいのは、Bが死刑執行されると告げられた後のAが恩赦を受ける確率P(A=o|B=t)である。
ベイズの公式により
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o)
P(B=t|B=o)=0 Bが恩赦を受けるときBが死刑執行されると告げられる確率=0
P(B=t|C=o)=1 Cが恩赦を受けるときBが死刑執行されると告げられる確率=1
問題は P(B=t|A=o)
恩赦を受けるのがAであるときに看守がCではなくBが死刑執行されると告げる確率は示されていない。
この確率をｐとすると
P(A=o|B=t)は p/(p+1)となる。
もちろんp=0.5であれば、P(A=o|B=t)=1/3と看守に告げられる前と同じである。
ここでｐが一様分布からさまざなβ分布に従うとするとどうなるか、グラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
左の緑が看守がBとCが死刑執行予定であるときにBを選んで答える確率分布。
右の青が看守がBと告げたときのAが恩赦を受ける確率の分布。
370:卵の名無しさん
18/01/03 21:28:54.72 CFofwzsi.net
無情報分布として一様分布を考えると
Aが恩赦を受ける確率の期待値（平均値）は
> 1-log(2)
[1] 0.3068528
となる。
p/(p+1)を [0,1]で定積分すれば求まる。
無情報分布として一様分布を考えると
Aが恩赦を受ける確率の期待値（平均値）は
> 1-log(2)
[1] 0.3068528
となる。
Cが恩赦を受ける確率の期待値（平均値）は
> log(2)
[1] 0.6931472
当然、Bが恩赦を受ける確率は０ 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)
371:卵の名無しさん
18/01/04 16:17:58.24 nNHCcGvL.net
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o) = (p + q+ 1-q)*1/3 = (p + 1)*1/3
P(B=t|A=o)*P(A=o) = p*1/3
P(B=t|B=o)*P(B=o) = q*1/3
P(B=t|C=o)*P(C=o) = (1-q)*1/3
P(A=o|B=t) = p / (p +1)
372:卵の名無しさん
18/01/05 09:01:44.41 HIwtxZA8.net
> vonNeumann
function(PDF,xmin=0,xmax=1,N=10000,Print=TRUE,...){
xx=seq(xmin,xmax,length=N+1)
ymax=max(PDF(xx))
Ux=runif(N,xmin,xmax)
Uy=runif(N,0,ymax)
Rand=Ux[which(Uy<=PDF(Ux))]
if(Print){
hist(Rand,xlim=c(xmin,xmax),freq=FALSE,col=sample(colors(),1),main='',...)
AUC=integrate(PDF,xmin,xmax)$value
lines(xx,sapply(xx,function(x)PDF(x)/AUC))
}
hdi=HDInterval::hdi(Rand)
print(c(hdi[1],hdi[2]),digits=4)
invisible(Rand)
}
373:卵の名無しさん
18/01/05 13:47:38.14 2i23lFVy.net
####
stanString2='
functions{
real jisaku_lpdf(real y, real a, real b){
real temp;
temp = a*y + b;
return log(temp);
}
}
data{
real a;
real b;
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real q;
q = p/(p+1);
}
model{
target += jisaku_lpdf( p | a, b);
}
'
374:卵の名無しさん
18/01/06 05:46:40.68 Qus13FVj.net
amnesty <- function(p) p/(p+1)
f=function(ab,k=10^4,ma=1/3){
a=ab[1]
b=ab[2]
rn=rbeta(k,a,b)
(mean(amnesty(rn)) - ma)^2
}
optim (c(1,1),f, method='Nelder')
optim (c(2,2
375:),f, method='Nelder')
376:卵の名無しさん
18/01/07 20:56:03.58 63fon0Wx.net
■3囚人問題（英: Three Prisoners problem）
ある監獄にA、B、Cという3人の囚人がいます
3人のうちランダムに選ばれた1人に恩赦が出ます
誰が死刑になるかは看守は本人には答えない
囚人Aに看守が「Bは死刑になる」と教えた
この時、看守は一定の確率で嘘をつく。
BとCが死刑になるときは50%の確率でBと答える。
囚人Aに恩赦が与えられる確率は何％でしょうか？
377:卵の名無しさん
18/01/09 19:32:17.33 WXFR1iXr.net
■3囚人問題（英: Three Prisoners problem）
ある監獄にA、B、Cという3人の囚人がいます。
3人のうちランダムに選ばれた1人に恩赦が出ます。
誰が死刑になるかは看守は決して本人には教えない。
囚人AがB、Cの少なくともどちらかは死刑になるのだから教えてくれと看守に尋ねた。看守が本人に教えるのではないので「Bは死刑になる」とAに教えた。
この時、看守は1/3の確率で嘘をつく。
BとCが死刑になるときは50%の確率でBと答える。
囚人A、B、Cに恩赦が与えられる確率はそれぞれ何％でしょう？
378:卵の名無しさん
18/01/09 19:41:09.46 WXFR1iXr.net
>>447
p : P(B=t|A=o)Aが恩赦（BとCが死刑執行される）とき看守がBと答える確率
q : 看守が嘘をつく確率
P(B=t|B=o) Bが恩赦を受けるときBが死刑執行されると告げられる確率 = q
P(B=t|C=o) Cが恩赦を受けるときBが死刑執行されると告げられる確率 = 1-q
P(A=o|B=t) = P(B=t|A=o)*P(A=o) / P(B=t)
P(B=t) = P(B=t|A=o)*P(A=o) + P(B=t|B=o)*P(B=o) + P(B=t|C=o)*P(C=o)
= p * P(A=o) + q * P(B=o) + (1-q) * P(C=o)
P(A=o|B=t) = p*P(A=o) / ( p*P(A=o) + q * P(B=o) + (1-q) * P(C=o) )
P(A=o)= P(B=o)= P(C=o) = 1/3ならば
P(A=o|B=t) = p /(p+1)
P(B=o|B=t) = q/(p+1)
P(C=o|B=t) = (1-q)/(p+1)
379:卵の名無しさん
18/01/18 15:55:28.61 mFfO4JsF.net
プレーヤーが選んだ箱をA、司会者モンティーホールが開けたハズレの箱をB、残った箱をCとする。
Aがアタリ(atari)の確率をP(A=a)
司会者がBを開ける(open)確率をP(B=o)と表すことにする。
残った箱Cがアタリである確率P(C=a|B=o)は
ベイズの公式から
P(C=a|B=o) = (P(B=o|C=a)P(C=a)) / P(B=o)
P(B=o) = P(B=o|A=a)P(A=a) + P(B=o|B=a)P(B=a) + P(B=o|C=a)P(C=a)
= P(B=o|A=a)*1/3 + 0*1/3 + 1*1/3 ここで P(B=o|A=a)=pとおくと
= p*1/3 + 1/3
ゆえに
P(C=a|B=o) = (P(B=o|C=a)P(c=a)) / P(B=o) = (1*1/3) / (p*1/3 + 1/3) = 1/(p+1)
380:卵の名無しさん
18/01/18 15:58:49.78 mFfO4JsF.net
BがハズレとわかったあとでAがアタリである確率
P(A=a|B=o) = P(B=o|A=a)P(A=a)/P(B=o)
P(B=o) = P(B=o|A=a)P(A=a) + P(B=o|B=a)P(B=a) + P(B=o|C=a)P(C=a)
P(B=o|A=a)はAがアタリであるときにBがハズレとして開けられる確率pは問題で示されていない。
不十分理由の原則に準じてpを０．５とするか一様分布に従うとするのが一般的だと思う。
P(B=o|A=a)=pとおくと
P(B=o) = P(B=o|A=a)P(A=a) + P(B=o|B=a)P(B=a) + P(B=o|C=a)P(C=a)
= p*1/3 + 0*1/3 + 1*1/3
= p*1/3 +1/3
ゆえに
P(A=a|B=o) = (p*1/3) / ( p*1/3 + 1/3 ) = p/(p+1)となる。
p=0.5ならBがハズレというデータはAがあたりの確率に影響を与えず1/3である。
381:卵の名無しさん
18/01/18 15:59:07.36 mFfO4JsF.net
p: Aがアタリの時に司会者がBを開ける確率
P(A=a|B=o) = p/(p+1) Bが開けられた後、Aがアタリの確率 (1)
P(C=a|B=o) = 1/(p+1) Bが開けられた後、Cがアタリの確率 (2)
(1)/(2) = p なので (2)は(1)以上である。(∵0<= p <=1)
ゆえに
残った箱Cの方がアタリの確率は高い。
382:卵の名無しさん
18/01/21 07:14:29.94 VjwevIsc.net
ゴルゴ15は1発1中
とする。
各々10000発撃ったときゴルゴ15の命中数の期待値はいくらか？
確率密度とかベータ分布とかを使わずに説明するなら、重み付き平均という考え方で説明するしかないかな？
命中率が0.5なら2回に1回は1発1中（確率0.5)
命中率が0.8なら10回に8回は1発1中（確率0.8)
となる。
体重100ｋｇの牛が100頭
体重99kgの牛が99頭
体重98kgの牛が98頭
・・・
体重2kgの牛が2頭
体重1kgの牛が1頭
牛の平均体重の計算と同じ
n=100
x=seq(0,1,length=n+1)
sum(x*x/sum(x))
sum(x^2)/sum(x)
2/3
(sum_x=n*(n+1)/2/n) # (n+1)/2
(sum_x2=n*(n+1)*(2*n+1)/6/(n^2)) # (n+1)*(2*n+1)/n/6
sum_x2/sum_x # (2*n+1)/n/3 = 2/3+1/3/n
n→∞
で2/3に集束する。 命中数の期待値は10000*2/3=6667
383:卵の名無しさん
18/03/08 13:29:28.47 1C4m5UpT.net
URLﾘﾝｸ(i.imgur.com)
リリースされているデータだと効き目はこんなもん
384:卵の名無しさん
18/03/08 17:51:23.53 jJMnC74H.net
R でログランク検定を実行するには，観察時間を示す変数を time，打ち�
385:ﾘりフラグを event，グループを group とすれば，survdiff(Surv(time,event)~group) とすればよい。この例の場合なら，下枠内の通り。 なお，一般化ウィルコクソン検定をするには survdiff(Surv(time,event)~group,rho=1) とすればよい。 it13-4-2006.R require(survival) time <- c(4,6,8,9,5,7,12,14) event <- c(1,1,1,1,1,1,1,1) group <- c(1,1,1,1,2,2,2,2) survdiff(Surv(time,event)~group)
386:卵の名無しさん
18/03/08 21:28:51.57 jJMnC74H.net
サイコロを３０００回振って１の目が４９０回でたサイコロはイカサマサイコロか？
１０％までの歪みは許容する。
387:卵の名無しさん
18/03/08 21:29:24.44 jJMnC74H.net
crooked <- function(n,r,H0=0.5,d=0.1,credMass=0.95,Print=TRUE,...){
hdi=HDInterval::hdi(qbeta,credMass,shape1=1+r,shape2=1+n-r)
return(hdi)
}
crooked(3000,490,H0=1/6,xlim=c(0.12,0.20))
binom::binom.confint(490,3000)
1/6*c(0.9,1.1)
グラフ化するとHDI ⊂ ROPEになっている。
URLﾘﾝｸ(i.imgur.com)
388:卵の名無しさん
18/03/12 07:04:43.41 RtLIiZVI.net
# URLﾘﾝｸ(www.ncbi.nlm.nih.gov)
test\disease present absent
pos TP(a) FN(b)
neg FP(c) TN(d)
a+c=33
b+d=15
a+b=33
a/(a+c)=0.697
d/(b+d)=0.333
a/33=0.697
a=23
d/15=0.333
d=5
c=10
b=10
test\disease present absent
pos 23 TP(a) 10 FN(b)
neg 10 FP(c) 5 TN(d)
# URLﾘﾝｸ(statpages.info)
389:卵の名無しさん
18/03/12 13:24:27.62 UNJR7sdw.net
頭がある という所見は髄膜炎の診断に感度１００％である。
角がある という所見は髄膜炎の診断に特異度１００％である。
俺の経験上、こういう話をシリツ医大卒に振っても興味を示す奴はいないね。
390:卵の名無しさん
18/03/20 23:04:27.14 ql3rAyFF.net
pLR=.99/.01
nLR=.01/.99
preo=0.001/(1-0.001)
preo=0.1/(1-0.1)
poso=preo*pLR
poso/(1+poso)
391:卵の名無しさん
18/03/21 20:43:38.90 5ZNg8kVw.net
# modified Wald by Agresti and Coull
mWald <- function(S,n,cl=0.95){
z=qnorm(1-(1-cl)/2)
n_tilde=n+z^2
p_tilde=(S+z^2/2)/n_tilde
W=z*sqrt(p_tilde*(1-p_tilde)/n_tilde)
CI=c(mean=S/n,lower=p_tilde-W,upper=p_tilde+W)
return(CI)
}
binom::binom.agresti.coull(S,n)
392:卵の名無しさん
18/03/22 08:58:41.40 IaVjTmR+.net
3/n ≒ 1 - n√0.05
393:卵の名無しさん
18/03/23 10:37:59.90 OxTHnFDc.net
#標準偏差の信頼区間
sdCI=function(x,conf.level=0.95){
n=length(x)
df=n-1
sd=sd(x)
lwr=qchisq((1-conf.level)/2,df)
upr=qchisq((1-conf.level)/2,df,lower.tail=FALSE)
CI=data.frame(lower=sqrt((n-1)/upr)*sd,sd=sd,upper=sqrt((n-1)/lwr)*sd )
return(CI)
}
##差の信頼区間(生データなし)
DifCI=function(n1,n2,m1,m2,sd1,sd2){
pooledV=((n1-1)*sd1^2+(n2-1)*sd2^2)/(n1-1+n2-1)
SE12=sqrt((1/n1+1/n2)*pooledV)
w=qt(.975,n1-1+n2-1)*SE12
ci=c(m1-m2-w,m1-m2+w)
names(ci)=c("lower","upper")
return(ci)
}
# t検定(生データなし)
T.test=function(n1,n2,m1,m2,sd1,sd2){
SE12=sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
T=(m1-m2)/SE12
pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
394:卵の名無しさん
18/03/23 10:38:50.82 OxTHnFDc.net
## m sd SEM n
non=c(0.52,0.25,0.027,88)
rec=c(0.38,0.32,0.034,89)
eli=c(0.40,0.26,0.048,28)
# 生データなしで分散分析
lh=rbind(non,rec,eli)
colnames(lh)=c("m","sd","SEM","n") ; lh
mean.G=sum(lh[,"m"]*lh[,"n"])/sum(lh[,"n"])
SS.bit=sum((lh[,"m"]-mean.G)^2*lh[,"n"])
SS.wit=sum(lh[,"sd"]^2*(lh[,"n"]-1))
df.bit=nrow(lh)-1
df.wit=sum(lh[,"n"]-1)
MS.bit=SS.bit/df.bit
MS.wit=SS.wit/df.wit
F.ratio=MS.bit/MS.wit
pf(F.ratio,df.bit,df.wit,lower.tail=FALSE) # 0.003720507
(η2=(SS.bit)/(SS.bit+SS.wit)) # 0.05387927
395:卵の名無しさん
18/03/23 10:40:13.81 OxTHnFDc.net
#prop.test(S,n)$conf.int # S:success, n:number of trial
prop=function(S,n)prop.test(S,n)$conf.int[1:2]
#binom.test(S,n)$conf.int # S:success, n:number of trial
binom=function(S,n) binom.test(S,n)$conf.int[1:2]
## S:success, n:number of trial
binomCI=function(S,n,conf.level=0.95){
upper=uniroot(f=function(x)pbinom(S,n,x)-(1-conf.level)/2,c(0,1))$root
lower=uniroot(f=function(x)pbinom(S,n,x,lower.tail=FALSE)-(1-conf.level)/2,c(0,1))$root
CI=data.frame(lower=lower, upper=upper)
return(round(CI,3))}
##
wCI=function(S,n,cl=0.95){ #modified Wald CI for Proportion
z=qnorm(1-(1-cl)/2)
p.=(S+0.5*z^2)/(n+z^2)
W=z*sqrt((p.*(1-p.))/(n+z^2))
C.I.=data.frame(lower=p.-W,upper=p.+W)
return(C.I.)
}
##
bCI<-function(r,n,cl=0.95){
p<-r/n
q<-1-p
sd<-sqrt (n*p*q)
SE<-sd/sqrt(n)
Z<-qnorm(1-(1-cl)/2)
ci<-data.frame(lower=p-Z*SE/sqrt(n),upper=p+Z*SE/sqrt(n))
return(ci)
}
396:卵の名無しさん
18/03/25 14:44:59.55 Z8Ybstvj.net
Signif?cosis. Manifested by a failure to discern between
biological and statistical signif?cance (6). Individuals with
signif?cosis fail to realize that just because something is
unlikely to have occurred by chance doesn’t mean it’s
important (17). See also Borderline Probability Disorder.
Borderline Probability Disorder. Afflicted individuals
may dismiss the potential importance of results with P=0.06
while unquestioningly accepting the importance of results with P=0.05
URLﾘﾝｸ(journals.asm.org)
397:卵の名無しさん
18/03/25 16:56:05.47 Z8Ybstvj.net
r1=10 ; r2=16
n1=5936 ; n2=3403
a=r1
b=n1-r1
c=r2
d=n2-r2
Fragile <- function(a,b,c,d){
mat=matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list(c("Treated","Control"),c("Event","No Event")))
print(mat)
p=fisher.test(mat)$p.value
cat(paste('p-value =',round(p,5)))
i=0
while(p < 0.05 ){
p=fisher.test(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE))$p.value
i=i+1
}
cat('\n','\n')
cat(paste('Fragile Index = ',i,'\n','\n'))
print(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE,dimnames=list(c("Treated","Control"),c("Event","No Event"))))
cat(paste('p-value =',round(p,5)))
invisible(i)
}
398:卵の名無しさん
18/03/26 07:18:15.20 HyN4R5wr.net
# High Normal Total
# Event a b x
# No Event c d y
sqrt(1.37*4.09)
x+y=43407
x=1280
y=43407-1280
HR=x/(x+y) ; HR
HR1= 2.37*HR ; HR1
a/(a+c)=HR1
c=(1-HR1)*a/HR1
V(logRR)=1/a-1/x+1/c-1/y
SE(logRR)=sqrt(1/a-1/x+1/c-1/y)
HR1*exp(1.96*SE(logRR))=4.09
399:卵の名無しさん
18/03/26 07:19:05.61 HyN4R5wr.net
library(fmsb)
hazardratio <-function (a, b, PT1, PT0, conf.level = 0.95) {
.M <- a + b
.T <- PT1 + PT0
.MAT <- matrix(c(a, b, .M, PT1, PT0, .T), 3, 2)
colnames(.MAT) <- c("Cases", "Person-time")
rownames(.MAT) <- c("Exposed", "Unexposed", "Total")
class(.MAT) <- "table"
# print(.MAT)
ESTIMATE <- (a/PT1)/(b/PT0)
norm.pp <- qnorm(1 - (1 - conf.level)/2)
.CHI <- (a - (PT1/.T) * .M)/sqrt(.M * (PT1/.T) * (PT0/.T))
p.v <- 2 * (1 - pnorm(abs(.CHI)))
RRL <- ESTIMATE * exp(-norm.pp * sqrt(1/a + 1/b))
RRU <- ESTIMATE * exp(norm.pp * sqrt(1/a + 1/b))
CINT <- c(RRL, RRU)
attr(CINT, "conf.level") <- conf.level
RVAL <- list(p.value = p.v, conf.int = CINT, estimate = ESTIMATE,
method = "Incidence rate ratio estimate and its significance probability",
data.name = paste(deparse(substitute(a)), deparse(substitute(b)),
deparse(substitute(PT1)), deparse(substitute(PT0))))
class(RVAL) <- "htest"
return(RVAL)
}
f <- function(a,u0=1.37) hazardratio(a, x, a/HR1, x+y)$conf[1] - u0
uniroot(f,c(1,1000))$root
a=13
(1-HR1)*a/HR1
c=173
400:卵の名無しさん
18/03/26 16:51:22.07 xQgYtict.net
Run any R code you like. There are over three thousand R packages preloaded.
URLﾘﾝｸ(rdrr.io)
401:卵の名無しさん
18/03/27 06:27:14.34 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
Y=rnorm(N,1,1)
n=17
f <- function(){
x=sample(X,n)
y=sample(Y,n)
x-y
}
d=replicate (10^5, f())
hist(d)
g <- function(){
x=sample(X,n)
y=sample(Y,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
mean(p<0.05)
402:卵の名無しさん
18/03/27 06:36:24.73 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
Y=rnorm(N,1,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(Y,n)
x-y
}
d=replicate (10^5, f())
hist(d)
g <- function(){
x=sample(X,n)
y=sample(Y,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))
403:卵の名無しさん
18/03/27 06:42:13.61 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
Y=rnorm(N,1,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(Y,n)
mean(x)-mean(y)
}
d=replicate (10^5, f())
hist(d)
sd(d)
g <- function(){
x=sample(X,n)
y=sample(Y,n)
t.test(x,y)$p.value
}
p=replicate (10^5, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))
404:卵の名無しさん
18/03/27 06:42:36.33 7KQq5UQg.net
URLﾘﾝｸ(rsos.royalsocietypublishing.org)
405:卵の名無しさん
18/03/27 06:50:47.55 7KQq5UQg.net
> N=1000
> X=rnorm(N,0,1)
> Y=rnorm(N,1,1)
> n=16
> f <- function(){
+ x=sample(X,n)
+ y=sample(Y,n)
+ mean(x)-mean(y)
+ }
> d=replicate (10^5, f())
> hist(d)
> sd(d)
[1] 0.3574393
>
> g <- function(){
+ x=sample(X,n)
+ y=sample(Y,n)
+ t.test(x,y)$p.value
+ }
> p=replicate (10^5, g())
> hist(p)
> hist(-log10(p))
> mean(log10(p)<log10(0.05))
[1] 0.79713
>
>
406:卵の名無しさん
18/03/27 06:52:28.91 7KQq5UQg.net
> power.t.test(n=16,delta=1)
Two-sample t test power calculation
n = 16
delta = 1
sd = 1
sig.level = 0.05
power = 0.7813965
alternative = two.sided
NOTE: n is number in *each* group
407:卵の名無しさん
18/03/27 08:07:28.48 7KQq5UQg.net
N=1000
X=rnorm(N,0,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(X,n)
mean(x)-mean(y)
}
d=replicate (10^5, f())
hist(d)
sd(d)
g <- function(){
x=sample(X,n)
y=sample(X,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))
408:卵の名無しさん
18/03/27 10:32:56.09 7KQq5UQg.net
>>385
理論上、
mean(p<0.05)は
power.t.test(n=16,delta=1)$powerに一致。
409:卵の名無しさん
18/03/28 09:07:02.26 lnrrtf8h.net
calc.FPR = function(nsamp,pval,sigma,prior,delta1){
# sdiff=sqrt(sigma^2/nsamp + sigma^2/nsamp)
ns1=nsamp
ns2=nsamp
# CALC sdiff = sd of difference between means
sdiff=sqrt(sigma^2/ns1 + sigma^2/ns2)
df=ns1+ns2-2
# Note FPR doesn't need calculation of power for p-equals case
#
#under H0, use central t distribution
tcrit=qt((1-pval/2),df,ncp=0)
x0=tcrit
y0=dt(x0,df,0)
#
# under H1 use non-central t distribution
ncp1=delta1/sdiff #non-centrality paramater
x1=x0 #tcrit
y1=dt(x1,df,ncp=ncp1)
# check solution
# pcheck=pt(y1,df,ncp=ncp1)
# pcheck
# Calc false positive risk
p0=2*y0
p1=y1
FPR=((1-prior)*p0)/(((1-prior)*p0) + prior*p1)
FPR
output=c(FPR,x0,y0,x1,y1)
return(output)
}
# end of function calc.FPR
#
410:卵の名無しさん
18/03/28 09:07:17.56 lnrrtf8h.net
# calc.FPR0 gives FPR for given nsamp, for p-less-than case
calc.FPR0 = function(nsamp,pval,sigma,prior,delta1){
myp=power.t.test(n=nsamp,sd=sigma,delta=delta1,sig.level=pval,type="two.sample",alternative="two.sided",power=NULL)
power = myp$power
PH1=prior
PH0=1-PH1
FPR0=(
411:pval*PH0/(pval*PH0 + PH1*power)) output=c(FPR0,power) return(output) } #
412:卵の名無しさん
18/03/28 09:27:44.91 lnrrtf8h.net
# False Positive Report Probability
FPRP <- function(n, p.value, prior, effect.size=1){
power=power.t.test(n=n,delta=effect.size,sig.level=p.value)$power
FP=(1-prior)*p.value
TP=prior*power
return(FP/(TP+FP))
}
FPRP(16,p.value=0.045,prior=0.10)
# サンプルサイズが違うとき
FPRP2 <- function(n1,n2, p.value, prior, effect.size=1,...){
power=pwr::pwr.t2n.test(n1, n2, d=effect.size, sig.leve=p.value, ...)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
FPRP2(10,6,p.value=0.05,prior=0.05)
413:卵の名無しさん
18/03/28 10:44:32.63 lnrrtf8h.net
# False Positive Report Probability
FPRP <- function(n, p.value, prior, effect.size=1){
power=power.t.test(n=n,delta=effect.size,sig.level=p.value)$power
FP=(1-prior)*p.value
TP=prior*power
return(FP/(TP+FP))
}
FPRP(16,p.value=0.045,prior=0.10)
# サンプルサイズが違うとき
FPRP2 <- function(n1,n2, p.value, prior, effect.size=1,...){
power=pwr::pwr.t2n.test(n1, n2, d=effect.size, sig.leve=p.value, ...)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
FPRP2(10,6,p.value=0.05,prior=0.05)
414:卵の名無しさん
18/03/28 18:20:29.83 izUkh1ln.net
更新版の
URLﾘﾝｸ(rsos.royalsocietypublishing.org)
を読んでたらRのスクリプトが公開されていた。
URLﾘﾝｸ(figshare.com)
415:卵の名無しさん
18/03/28 20:36:10.95 izUkh1ln.net
FPRP = P(H0|Rejected) = P(Rejected|H0)P(H0) / [ P(Rejected|H0)P(H0) + P(Rejected|H1)P(H1)]
= alpha*(1-prior) / [ alpha*(1-prior) + power*prior ] (1)
= p.value*(1-prior) / [ p.value*(1-prior) + power*prior ] (2)
416:卵の名無しさん
18/03/29 14:32:00.14 HGAAnleD.net
シミュレーションデータ
> length(CONTROLS) ; length(HYPERTENSION)
[1] 17
[1] 18
> mean(CONTROLS) ; mean(HYPERTENSION)
[1] 263
[1] 257
> sd(CONTROLS) ; sd(HYPERTENSION)
[1] 87
[1] 59
で条件に合致。
> t.test(CONTROLS,HYPERTENSION,var=TRUE)
Two Sample t-test
data: CONTROLS and HYPERTENSION
t = 0.24003, df = 33, p-value = 0.8118
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-44.85718 56.85718
JAGSを使ってMCMCしてみた(パッケージBESTを改造して使用）。
URLﾘﾝｸ(i.imgur.com)
理論値より信頼区間幅が広くなるが、まあ、シミュレーションなのでこんなもんだろう。
> HDIofMCMC(muDiff)
[1] -48.86974 61.99241
自分で計算してグラフ化しながら読み進むのは楽しい。
417:卵の名無しさん
18/03/29 22:09:13.07 HGAAnleD.net
library(pwr)
args(pwr.t2n.test)
n1=17;n2=18;m1=263;m2=257;sd1=87;sd2=59
sp_sq = ((n1-1)*sd1^2 + (n2-1)*sd2^2)/(n1-1+n2-1) # pooled variance
sp=sqrt(sp_sq) ; sp # pooled sd
sep=sqrt(1/n1+1/n2)*sp # se for pooled sd
f20.1 <- function(x) pwr.t2n.test(n1=17,n2=18, d=x, sig.level = 0.05)$power
xx=0:150
plot(xx,sapply(xx/sp,f20.1),type='l',ylab='power',xlab='mean difference')
f20.1(0)
f20.1(1)
#
z1=383;z2=373;n1=7685;n2=7596
p1=z1/n1
p2=z2/n2
p1/p2
ES.h(p1,p2)
g20.1 <- function(x){
pwr.2p2n.test(h=ES.h(p1,x*p1),n1=n1,n2=n2,sig.level=0.05)$power
}
g20.1(1)
g20.1(0.7)
rr=seq(0,1,by=0.01)
plot(rr,sapply(rr,g20.1),type='l',xlab='relative risk ratio',ylab='power')
418:卵の名無しさん
18/03/30 06:00:08.04 oNlAm/3j.net
# False Positive Report Probability for Proportion
FPRPP <- function(r1,r2,n1,n2, prior){
p1=r1/n1
p2=r2/n2
p.value=prop.test(c(r1,r2),c(n-1,n2))$p.value
pwr::pwr.2p2n.test(h=ES.h(p1,x*p1),n1=n1,n2=n2,sig.level=p.value)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
r1=10 ; r2=16
n1=5936 ; n2=3403
FPRPP(10,16,5936,3403,0.5)
419:卵の名無しさん
18/03/30 06:14:57.20 oNlAm/3j.net
# False Positive Report Probability for Proportion
FPRPP <- function(r1,r2,n1,n2, prior){
p.value=prop.test(c(r1,r2),c(n1,n2))$p.value
p1=r1/n1
p2=r2/n2
h=pwr::ES.h(p1,p2)
power=pwr::pwr.2p2n.test(h=h,n1=n1,n2=n2,sig.level=p.value)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
return(FP/(TP+FP))
}
r1=10 ; r2=16
n1=5936 ; n2=3403
FPRPP(r1,r2,n1,n2, prior=0.5)
420:卵の名無しさん
18/03/30 07:49:47.13 oNlAm/3j.net
URLﾘﾝｸ(clincalc.com)
421:卵の名無しさん
18/03/30 09:44:46.63 oNlAm/3j.net
# Fragile Index
FragileIndex<- function(r1,r2,n1,n2){
if(
422:r1/n1 < r2/n2){ a=r1 b=n1-r1 c=r2 d=n2-r2 }else{ a=r2 b=n2-r2 c=r1 d=n1-r1 } mat=matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list(c("low","high"),c("Event","No Event"))) print(mat) p=fisher.test(mat)$p.value cat(paste('p-value =',round(p,5))) i=0 while(p < 0.05 ){ p=fisher.test(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE))$p.value i=i+1 } cat('\n','\n') cat(paste('Fragile Index = ',i,'\n','\n')) print(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE,dimnames=list(c("low+FI","high-FI"),c("Event","No Event")))) cat(paste('p-value =',round(p,5))) invisible(i) }
423:卵の名無しさん
18/03/30 10:30:37.99 oNlAm/3j.net
# URLﾘﾝｸ(www.jclinepi.com)(13)00466-6/pdf
# Fragile Index
FragileIndex<- function(r1,r2,n1,n2,Fisher=TRUE){
if(r1/n1 < r2/n2){
a=r1
b=n1-r1
c=r2
d=n2-r2
}else{
a=r2
b=n2-r2
c=r1
d=n1-r1
}
mat=matrix(c(a,b,c,d),2,2,byrow=TRUE,dimnames=list(c("low","high"),c("Event","No Event")))
print(mat)
FUN=ifelse(Fisher,fisher.test,chisq.test)
p=FUN(mat)$p.value
cat(paste('p-value =',round(p,5)))
i=0
while(p < 0.05 ){
p=fisher.test(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE))$p.value
i=i+1
}
cat('\n','\n')
cat(paste('Fragile Index = ',i,'\n','\n'))
print(matrix(c(a+i,b-i,c,d),2,2,byrow=TRUE,dimnames=list(c("low+FI","high-FI"),c("Event","No Event"))))
cat(paste('p-value =',round(p,5)))
invisible(i)
}
424:卵の名無しさん
18/03/30 14:34:22.51 oNlAm/3j.net
rm(list=ls(all=TRUE))
graphics.off()
source('URLﾘﾝｸ(raw.githubusercontent.com)')
BMV_vs_ETI <- function(zi=44 , Ni=1018 , zj=43 , Nj=1022, ROPE=c(-0.005,0.005),COLOR='skyblue'){
y=c(rep(1,zi),rep(0,Ni-zi),rep(1,zj),rep(0,Nj-zj))
s=rep(1:2,c(Ni,Nj))
shape1=1 ; shape2=1 # JAGS prior : beta(shape1,shape2)
myData=data.frame(y=y,s=s)
Ntotal = length(y)
Nsubj = length(unique(s))
dataList = list(
y = y ,
s = s ,
Ntotal = Ntotal ,
Nsubj = Nsubj
)
# JAGS model
modelString = paste0("
model {
for ( i in 1:Ntotal ) {
y[i] ~ dbern( theta[s[i]] )
}
for ( sIdx in 1:Nsubj ) {
theta[sIdx] ~ dbeta(", shape1,',' , shape2," )
}
}
")
# end of modelString
425:卵の名無しさん
18/03/30 14:34:54.83 oNlAm/3j.net
writeLines( modelString , con="TEMPmodel.txt" )
jagsModel = jags.model( file="TEMPmodel.txt" , data=dataList)
update( jagsModel)
codaSamples = coda.samples( jagsModel , variable="theta", n.iter=10000 )
mcmcMat=as.matrix(codaSamples)
dif=mcmcMat[,1]-mcmcMat[,2]
plotPost(dif,ROPE=ROPE,compVal=0,xlab=quote(Delta),breaks=30,cenT='mean',col=COLOR)
invisible(dif)
}
# favorable functional survival at day 28
BMV_vs_ETI(zi=44 , Ni=1018 , zj=43, Nj=1022, ROPE=c(-0.01,0))
# Survival to hospital admission 294/1018 vs 333/1022
BMV_vs_ETI(294,1018,333,1022,ROPE=c(-0.01,0.01))
# global survival at day 28 55/1018 vs 54/1022
BMV_vs_ETI(55,1018,54,1022,ROPE=c(-0.01,0.01))
426:卵の名無しさん
18/03/30 22:20:19.00 /bcqdZEY.net
Generic=rlnorm(10,4.5,0.2)
Brand=rlnorm(20,4.3,0.1)
while(!Equiv(Generic,Brand)){
Generic=rlnorm(10,4.5,0.2)
Brand=rlnorm(20,4.3,0.1)
}
Generic
Brand
t.test(Generic,Brand)$p.value
tt=t.test(Generic,conf.level = 0.90) ; tt
lo=tt$conf[1]
hi=tt$conf[2]
m=mean(Brand)
0.80*m < lo
hi < 1.25*m
427:卵の名無しさん
18/03/31 07:28:51.22 7NaQ5iyd.net
# How big the sample should be to show non-inferiority
r1=44; n1=1018;r2=43;n2=1022
BMV_vs_ETI(r1,n1,r2,n2)
risk.difference(r1,r2,n1,n2)
f <- function(x,margin=0.01){
rd=fmsb::riskdifference(x*r1,x*r2,x*n1,x*n2)
lo=rd$conf[1]
hi=rd$conf[2]
- margin - lo
}
m=uniroot(f,c(0.001,5))$root ; m
BMV_vs_ETI(round(m*r1),round(m*n1),round(m*r2),round(m*n2),COLOR = 'pink')
428:卵の名無しさん
18/04/01 20:44:41.26 yXplTzVs.net
Freedman50 <- function(){
X=matrix(rnorm(50*100),ncol=100) # X[50 variables,100 samples]
Y=rnorm(100)
ans=lm(Y~ X[1,]+X[2,]+X[3,]+X[4,]+X[5,]+X[6,]+X[7,]+X[8,]+X[9,]+X[10,]+X[11,]+X[12,]+X[13,]+X[14,]+X[15,]+X[16,]+X[17,]+X[18,]+X[19,]+X[20,]+
X[21,]+X[22,]+X[23,]+X[24,]+X[25,]+X[26,]+X[27,]+X[28,]+X[29,]+X[30,]+X[31,]+X[32,]+X[33,]+X[34,]+X[35,]+X[36,]+X[37,]+X[38,]+X[39,]+X[40,]+
X[41,]+X[42,]+X[43,]+X[44,]+X[45,]+X[46,]+X[47,]+X[48,]+X[49,]+X[50,])
anova.ans=anova(ans)
p.value=anova.ans$`Pr(>F)`[1] # p.value
429: of regression smry=summary(ans) r.squared=smry$r.squared pv.coef=smry$coef[-1,4] pv.coef025=sum(pv.coef<0.25) # p.value of coefficient pv.coef005=sum(pv.coef<0.05) data.frame(p.value, r.squared,pv.coef025,pv.coef005) } Freedman50() res50=replicate(10^3,Freedman50()$p.value) BEST::plotPost(res50,compVal = 0.05)
430:卵の名無しさん
18/04/02 09:48:50.18 qPQfofGo.net
f<- function(d,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f(5/10)
power.t.test(delta=0.5,power=0.80)$n
dd=seq(0.01,0.99,by=0.01)
g<- function(x) f(x,.50)/f(x,.80)
sapply(dd,g)
(2*(abs(1.96+abs(qnorm(1-0.80))))^2)/(2*(abs(1.96+abs(qnorm(1-0.50))))^2)
(1.96+0.84)^2/(1.96+0)^2
431:卵の名無しさん
18/04/02 10:00:29.14 qPQfofGo.net
f<- function(d,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f(5/10)
power.t.test(delta=0.5,power=0.80)$n
dd=seq(0.01,0.99,by=0.01)
g<- function(x) f(x,.50)/f(x,.80)
sapply(dd,g)
(2*(abs(1.96+abs(qnorm(1-0.80))))^2)/(2*(abs(1.96+abs(qnorm(1-0.50))))^2)
(1.96+0.84)^2/(1.96+0)^2
dd=seq(0.01,0.99,by=0.01)
h<-function (x){
power.t.test(d=x, power=0.80)$n/power.t.test(d=x, power=0.5)$n
}
res=sapply (dd,h)
plot(dd,res)
summary (res)
432:卵の名無しさん
18/04/03 00:15:01.95 wYZjtgN6.net
FPRPP.F <- function(r1,r2,n1,n2, prior=0.5){
p.value=fisher.test(matrix(c(r1,n1-r1,r2,n2-r2),2))$p.value
p1=r1/n1
p2=r2/n2
power=statmod::power.fisher.test(p1,p2,n1,n2,alpha=p.value)
FP=(1-prior)*p.value
TP=prior*power
c(p.value=p.value,'false positive rate' = FP/(TP+FP))
}
r1=10;n1=19;r2=11;n2=12
FPRPP.F(r1,r2,n1,n2,0.5)
433:卵の名無しさん
18/04/03 14:26:09.59 wYZjtgN6.net
標準偏差の信頼区間
f247 <- function(n,conf.level=0.95){
df=n-1
lwr=qchisq((1-conf.level)/2,df)
upr=qchisq((1-conf.level)/2,df,lower.tail=FALSE)
c('Lower limit'=sqrt((n-1)/upr),'Upper limit'=sqrt((n-1)/lwr))
}
sd.ci=sapply(c(2,3,4,5,10,25,50,100,500,1000),f247)
colnames(sd.ci)=paste('n =',c(2,3,4,5,10,25,50,100,500,1000))
t(round(sd.ci,2))
434:卵の名無しさん
18/04/04 22:30:14.42 rn9pqwFF.net
URLﾘﾝｸ(www.rdocumentation.org)
π=P(H0=TRUE)=P(association)は
π=P(HA=TRUE)=P(association)の間違いじゃないだろうか？
435:卵の名無しさん
18/04/05 14:14:03.67 xAGVH6hI.net
N=100
xy=rnorm(N)
index=1:N
f254 <- function(r1=5,K=75,PRINT=TRUE){
x0=sample(index,r1) ; y0=sample(index,r1)
PV=numeric()
PV[1]=t.test(xy[x0] , xy[y0])$p.value
x.index=x0 ; y.index=y0
for(i in 2:(K-r1+1)){
add.x=sample(index[-x.index],1)
x.index=append(x.index,add.x)
add.y=sample(index[-y.index],1)
y.index=append(y.index,add.y)
PV[i]=t.test(xy[x.index],xy[y.index])$p.value
i=i+1
}
if(PRINT){
plot(r1:K,PV,type='l',lwd=2, log='y', ylim=c(0.001,1),
xlab='n',ylab='p-value')
abline(h=0.05,lty=3)
}
FPR=mean(PV<0.05)
return(FPR)
}
x11() ; par(mfrow=c(1,1)) ;f254()
x11() ;par(mfrow=c(3,3)) ; replicate(9,f254())
436:卵の名無しさん
18/04/05 14:43:46.32 xAGVH6hI.net
The simulated data were sampled from populations with Gaussian distribution and identical means and standard deviations, An unpaired t test was computed with n=5 in each group,
and the resulting P value is plotted at the left of the graph.Then the t test was repeated with one more value added to each group.
Then one more value was added to each group(n=7),and the P value was computed again. This was continued up to n=75.
こういうのを読むと自分でシミュレーションせずにはいられなくなる。
URLﾘﾝｸ(i.imgur.com)
Rのコードはこれ。
ｽﾚﾘﾝｸ(hosp板:410番)
437:卵の名無しさん
18/04/06 19:26:15.51 I1oriqts.net
Appendix Table 1. Joint probability of significance of test and truth of hypothesis
Truth of alternative hypothesis
Significance of test
Significant Not significant Total
True association (1 ) [True positive] [False negative]
No association -
(1 ) [False positive] (1 -
) (1 ) [True negative] 1
Total (1 )-
(1 ) (1 -
) (1 ) 1
440 COMMENTARY Journal of the National Cancer Institute, Vol. 96, No. 6, March 17, 2004
438:卵の名無しさん
18/04/07 13:38:05.89 4kW+Nj56.net
f<- function(d=1,beta=.80,alpha=.05) 2*(abs((qnorm(1-alpha/2))+abs(qnorm(1-beta)))/d)^2
f(5/10)
power.t.test(delta=0.5,power=0.80)$n
# power 50%->80%
dd=seq(0.01,0.99,by=0.01)
g<- function(x) f(x,.50)/f(x,.80)
plot(dd,sapply(dd,g))
(2*(abs(1.96+abs(qnorm(1-0.80))))^2)/(2*(abs(1.96+abs(qnorm(1-0.50))))^2)
(1.96+0.84)^2/(1.96+0)^2
dd=seq(0.01,0.99,by=0.01)
h<-function (x,pow1=0.80,pow2=0.50){
power.t.test(d=x, power=pow1)$n/power.t.test(d=x, power=pow2)$n
}
res=sapply (dd,h)
plot(dd,res)
summary (res)
439:卵の名無しさん
18/04/07 13:41:51.63 4kW+Nj56.net
N=10^5
D=rbeta(N,2,1)
K=rbeta(N,1,2)
par (mfrow=c(2,2))
hist(D)
hist(D-K)
hist (K)
hist(log(D/K))
quantile (D-K, probs=c(0.025,0.5,0.975))
quantile (D/K, probs=c(0.025,0.5,0.975))
summary (D-K)
summary (D/K)
exp(mean (log(D/K)))
median (D/K)
BEST::plotPost (D-K)
BEST::plotPost (D/K)
BEST::plotPost(log(D/K))
440:卵の名無しさん
18/04/07 13:46:36.71 4kW+Nj56.net
N=10^5
D=rbeta(N,2,1)
K=rbeta(N,1,2)
par (mfrow=c(2,2))
hist(D)
hist(D-K)
hist (K)
hist(log(D/K))
quantile (D-K, probs=c(0.025,0.5,0.975))
quantile (D/K, probs=c(0.025,0.5,0.975))
summary (D-K)
summary (D/K)
exp(mean (log(D/K)))
median (D/K)
NNT=1/(D-K)
BEST::plotPost (D-K)
BEST:: plotPost (NNT)
BEST::plotPost (D/K)
BEST::plotPost(log(D/K))
441:卵の名無しさん
18/04/07 14:23:06.29 4kW+Nj56.net
f3.6 <- function(delta,n){
pnorm(-qnorm(.975)-sqrt(n)*delta)+pnorm(qnorm(.975)-sqrt(n)*delta,lower=FALSE)
}
f3.6(0.6,9)
curve (f3.6(x,25),-2,2,xlab=quote(Delta), ylab='power')
curve (f3.6(x,16),add=TRUE)
curve (f3.6(x, 9),add=TRUE)
442:卵の名無しさん
18/04/08 11:09:58.60 vtyfqsbT.net
# HDI for chisq
curve(dchisq(x,9),0,30)
conf.level=0.95
n=10
f <- function(x,df=n-1,conf.level=0.95){
d <- function(d) pchisq(x+d,df) - pchisq(x,df) - conf.level
uniroot(d,c(0,df*3))$root
}
qchisq(1-conf.level,n-1)
low=seq(0,floor(qchisq(1-conf.level,n-1)),by=0.01)
plot(low,sapply(low,f),type='l',lwd=2)
opt=optimise(f,c(0,floor(qchisq(1-conf.level,n-1)))) ; opt
opt[[1]] ; opt[[1]]+opt[[2]]
pchisq(opt[[1]]+opt[[2]],n-1) - pchisq(opt[[1]],n-1)
qchisq(0.025,n-1) ; qchisq(0.975,n-1)
pchisq(qchisq(0.975,n-1),n-1)-pchisq(qchisq(0.025,n-1),n-1)
443:卵の名無しさん
18/04/08 13:30:22.84 vtyfqsbT.net
# 不偏分散u2が重要なのは（ランダムサンプリングでは）
# 不偏分散の期待値が母分散と一致するからです。
# 定理：E[u2]=σ2
N=1000
n=10
X=rpois(N3)
hist(X)
var(X)
k=10^4
mean(replicate(k,var(sample(X,n))))
VAR <- function(x){
n.x=length(x)
var(x)*(n.x-1)/n.x
}
mean(replicate(k,VAR(sample(X,n))))
444:卵の名無しさん
18/04/09 08:33:45.55 v1O7hrCx.net
# 6割以上正解 0.522
sum(dbinom(240:400,400,0.6))
binom.test(240,400,0.6, alt='greater')$p.value
# 6割未満正解 0.478
sum(dbinom(0:239,400,0.6))
binom.test(239,400,0.6,alt='less')$p.value
n=5
# 禁忌枝3個以上選択0.317
binom.test(3,n,1-0.6,alt='g')
sum(dbinom(3:n,n,1-0.6))
0.522*(1-0.317)
# 一般問題を1問1点、臨床実地問題を1問3点とし、
#（1）?（3）のすべての合格基準を満たした者を合格とする。
#（1）必修問題 160点以上/200点
#（2）一般問題及び臨床実地問題 208点以上/299点
#（3）禁忌肢問題選択数 3問以下
p=0.6
binom.test(160,200,p,alt='greater')$p.value
binom.test(200,290,p,alt='greater')$p.value
n=5
445:卵の名無しさん
18/04/09 11:00:00.53 8EAvIhhV.net
現状の国試の合格基準は
一般問題を1問1点、臨床実地問題を1問3点とし、
（1）～（3）のすべての合格基準を満たした者を合格とする。
（1）必修問題 160点以上/200点
（2）一般問題及び臨床実地問題 208点以上/299点
（3）禁忌肢問題選択数 3問以下
来年から400問に減るらしい。
合格基準の正解率は同じ、即ち、(160+208)/499*400=295以上正解を要し、かつ、禁忌肢選択も３問以下とする。
正解率80%の受験生集団の合格率が50%であった�
446:ﾆき 禁忌肢問題は何問出題されたか答えよ。 f<-function (n,p){ (a=binom.test(295,400,p,alt='greater')$p.value) (b=binom.test(n-3,n,p,alt='greater')$p.value) a*b } nn=3:40 p=0.80 plot (nn,sapply(nn, function (n)f(n,p)))
447:卵の名無しさん
18/04/09 20:23:30.66 wJm8c5Cc.net
# F ~ F(Φ,∞) => Φ*F ~ χ2(Φ1)
fp.100 <- function(){
N=10^4
phi=sample(1:100,1)
rF=rf(N,phi,Inf)
hist(rF,breaks=25,freq=FALSE,col=sample(colours(),1),main=paste('Φ1 = ',phi1,', Φ2 = Inf'),
xlab=' F ')
curve(df(x,phi,Inf),add=TRUE)
hist(phi*rF,breaks=25,freq=FALSE,col=sample(colours(),1),main=paste('Φ1 = ',phi1,', Φ2 = Inf'),
xlab=' Φ*F ')
curve(dchisq(x,phi),add=TRUE,lwd=2)
}
x11() ; par(mfrow=c(3,3)) ; for(i in 1:9) fp.100()
448:卵の名無しさん
18/04/10 08:46:18.28 sK+H3q9V.net
XOF <- function (shape, scale,n=375,med=53.5,lwr=48.0,upr=58.5){
f <- function (v){
sh=v[1]
sc=v[2]
xof=rgamma(n,sh,sc) # sh: shape, sc: scale
(median(xof)-med)^2+ (qgamma(0.025, sh,sc) - lwr)^2 + (qgamma(0.975,sh,sc) - upr)^2
}
opt=optim (c(shape,scale), f, method='N')$par
print (opt)
curve(dgamma(x,opt [1],opt[2]),0,100)
sim=rgamma (n,opt [1],opt[2])
print(quantile (sim, probs=c(0.025,0.5,0.975)))
invisible (opt)
}
XOF(400,7)
449:卵の名無しさん
18/04/10 10:02:47.35 sK+H3q9V.net
XOF <- function (shape, scale,n=375,med=53.5,lwr=48.0,upr=58.5){
f <- function (v){
sh=v[1]
sc=v[2]
N=10^4
xof=rgamma(N,sh,sc) # sh: shape, sc: scale
(median(xof)-med)^2+ (qgamma(0.025, sh,sc) - lwr)^2 + (qgamma(0.975,sh,sc) - upr)^2
}
opt=optim (c(shape,scale), f, method='N')$par
print (opt)
curve(dgamma(x,opt [1],opt[2]),0,100)
sim=rgamma (n,opt [1],opt[2])
print(quantile (sim, probs=c(0.025,0.5,0.975)))
invisible (opt)
}
xof=XOF(400,7,375,53.5,48.0,58.5)
tam=XOF(400,7,377,53.8,50.2,56.4)
X=rgamma (375,xof[1],xof[2])
T=rgamma (377,tam[1],tam[2])
wilcox.test (X,T)
t.test(X,T)
bss=10^3
f1 <- function (){
XO=sample (X,bss, replace=TRUE)
TA=sample (T,bss, replace=TRUE)
mean(XO-TA)
}
dif=replicate (10^3,f1())
quantile (dif,probs=c(.024,.5,.975))
URLﾘﾝｸ(imagizer.imageshack.com)
450:卵の名無しさん
18/04/10 12:35:44.49 sK+H3q9V.net
XOF <- function (shape, scale,n=375,med=53.5,lwr=48.0,upr=58.5){
f <- function (v){
sh=v[1]
sc=v[2]
N=10^4
xof=rgamma(N,sh,sc) # sh: shape, sc: scale
(median(xof)-med)^2+ (qgamma(0.025, sh,sc) - lwr)^2 + (qgamma(0.975,sh,sc) - upr)^2
}
opt=optim (c(shape,scale), f, method='L-BFGS-B')
opt_par=opt$par
print (opt_par)
if(opt$convergence!=0) print(opt)
x11()
curve(dgamma(x,opt_par[1],opt_par[2]),40,70,xlab='time (h)',ylab='',bty='l')
sim=rgamma (n,opt_par[1],opt_par[2])
print(quantile(sim, probs=c(.025,.5,.975)),digits=3)
invisible (opt_par)
}
xof=XOF(400,7,375,53.5,48.0,58.5) ; tam=XOF(1000,20,377,53.8,50.2,56.4)
X=rgamma (375,xof[1],xof[2]) ; T=rgamma (377,tam[1],tam[2])
wilcox.test (X,T)
t.test(X,T)
bss=10^3
f1 <- function (){
XO=sample (X,bss, replace=TRUE)
TA=sample (T,bss, replace=TRUE)
mean(XO-TA)
}
dif=replicate (10^3,f1())
print(quantile (dif,probs=c(.024,.5,.975)),digits=3)
451:卵の名無しさん
18/04/10 15:18:27.99 sK+H3q9V.net
xo=data.frame(n=375,med=53.5,l=48.0,u=58.5)
ta=data.frame(n=377,med=53.8,l=50.2,u=56.4)
(xo$l+xo$u)/2 ; (ta$l+ta$u)/2
sqrt(xo$l*xo$u) ; sqrt(ta$l*ta$u)
ci2sd.t <- function(n,L,U){
(U-L)/2/qt(.975,n-1)*sqrt(n)
}
ci2sd.n <- function(n,L,U){
(U-L)/2/qnorm(.975)*sqrt(n)
}
(sd.xo=ci2sd.n(xo$n,xo$l,xo$u)) ; ci2sd.t(xo$n,xo$l,xo$u)
(sd.ta=ci2sd.n(ta$n,ta$l,ta$u)) ; ci2sd.t(ta$n,ta$l,ta$u)
s1=sd.xo ; s2=sd.ta
n1=xo$n ; n2=ta$n
sd=sqrt(((n1-1)*s1^2 + (n2-1)*s2)/(n1 + n2 -2))
sd
# 両側検定
f7.11n <- function(n,delta,alpha=0.05){ # n1 = n2
df=n+n-2
ncp=delta/sqrt(1/n + 1/n) # ⊿= (μ1-μ0)/σ
pt(qt(alpha/2,df),df,ncp)+
pt(-qt(alpha/2,df),df,ncp,lower=FALSE)
}
f7.17n <- function(power,delta,alpha=0.05){
uniroot(function(n)f7.11n(n,delta,alpha)-power,c(3,10^6))$root
}
f7.17n(0.8,0.3/sd)
452:卵の名無しさん
18/04/10 15:27:03.10 sK+H3q9V.net
# 平均値信頼区間から必要なサンプルサイズを計算
xo=data.frame(n=375,med=53.5,l=48.0,u=58.5)
ta=data.frame(n=377,med=53.8,l=50.2,u=56.4)
(xo$l+xo$u)/2 ; (ta$l+ta$u)/2
sqrt(xo$l*xo$u) ; sqrt(ta$l*ta$u)
ci2sd.t <- function(n,L,U){
(U-L)/2/qt(.975,n-1)*sqrt(n)
}
ci2sd.n <- function(n,L,U){
(U-L)/2/qnorm(.975)*sqrt(n)
}
(sd.xo=ci2sd.t(xo$n,xo$l,xo$u)) ; ci2sd.n(xo$n,xo$l,xo$u)
(sd.ta=ci2sd.t(ta$n,ta$l,ta$u)) ; ci2sd.n(ta$n,ta$l,ta$u)
s1=sd.xo ; s2=sd.ta
n1=xo$n ; n2=ta$n
sd=sqrt(((n1-1)*s1^2 + (n2-1)*s2)/(n1 + n2 -2))
sd
453:卵の名無しさん
18/04/10 15:27:18.84 sK+H3q9V.net
# 両側検定
f7.11n <- function(n,delta,alpha=0.05){ # n1 = n2
df=n+n-2
ncp=delta/sqrt(1/n + 1/n) # ⊿= (μ1-μ0)/σ
pt(qt(alpha/2,df),df,ncp)+
pt(-qt(alpha/2,df),df,ncp,lower=FALSE)
}
f7.17n <- function(power,delta,alpha=0.05){
uniroot(function(n)f7.11n(n,delta,alpha)-power,c(3,10^6))$root
}
f7.17n(0.8,0.3/sd)
# 片側検定 H1: mean(xo) < mean(ta)
f7.16n <- function(n,delta,alpha=0.05){ # n1=n2
df=n+n-2
ncp=delta/sqrt(1/n+1/n) # ⊿= (μ1-μ0)/σ
pt(qt(alpha,df),df,ncp)
}
f7.19n <- function(power,delta,alpha=0.05){
uniroot(function(n)f7.16n(n,delta,alpha)-power,c(3,10^6))$root
}
f7.19n(0.80,-0.3/sd)
454:卵の名無しさん
18/04/10 20:33:17.78 sK+H3q9V.net
xo=c(455,med=53.7,49.5,58.5)
pl=c(230,med=80.2,72.6,87.1)
n1=xo[1] ;n2=pl[1]
ci2sd.t <- function(n,L,U){
(U-L)/2/qt(.975,n-1)*sqrt(n)
}
ci2sd.n <- function(n,L,U){
(U-L)/2/qnorm(.975)*sqrt(n)
}
pooled.sd <- function(sd1,sd2,n1,n2){
sqrt((sd1^2*(n1-1)+
455:sd2^2*(n2-1))/(n1+n2-2)) } sd1=ci2sd.t(n1,xo[3],xo[4]) ; sd1 sd2=ci2sd.t(n2,pl[3],pl[4]) ; sd2 sd.p=pooled.sd(sd1,sd2,n1,n2) ; sd.p pwr::pwr.t2n.test(n1,n2,d=-26.5/sd.p,alt='less') pwr::pwr.t2n.test(n1=NULL,n2=n2,d=-26.5/sd.p,power=0.80,alt='less') # t検定(生データなし,等分散不問,両側検定) Welch.test=function(n1,n2,m1,m2,sd1,sd2){ T=(m1-m2)/sqrt(sd1^2/n1+sd2^2/n2) df=(sd1^2/n1+sd2^2/n2)^2 / (sd1^4/n1^2/(n1-1)+sd2^4/n2^2/(n2-1)) p.value=2*pt(abs(T),df,lower.tail = FALSE) return(p.value) } m1=(xo[3]+xo[4])/2 ; m2=(pl[3]+pl[4])/2 Welch.test(n1,n2,m1,m2,sd1,sd2) T=(m1-m2)/sqrt(sd1^2/n1+sd2^2/n2) df=(sd1^2/n1+sd2^2/n2)^2 / (sd1^4/n1^2/(n1-1)+sd2^4/n2^2/(n2-1)) pt(T,df)
456:卵の名無しさん
18/04/11 11:22:50.41 KzjAFvro.net
# False Positive Report Probability for Proportion
FPRPP <- function(r1,r2,n1,n2, prior=0.5){
p.value=prop.test(c(r1,r2),c(n1,n2))$p.value
p1=r1/n1
p2=r2/n2
h=pwr::ES.h(p1,p2)
power=pwr::pwr.2p2n.test(h=h,n1=n1,n2=n2,sig.level=p.value)$power
FP=(1-prior)*p.value # p-less-than case
TP=prior*power
c(p.value=p.value,'false positive rate' = FP/(TP+FP))
}
457:卵の名無しさん
18/04/11 11:25:58.50 KzjAFvro.net
FPRPP.F <- function(r1,r2,n1,n2, prior=0.5){
p.value=fisher.test(matrix(c(r1,n1-r1,r2,n2-r2),2))$p.value
p1=r1/n1
p2=r2/n2
power=statmod::power.fisher.test(p1,p2,n1,n2,alpha=p.value)
FP=(1-prior)*p.value
TP=prior*power
c(p.value=p.value,'false positive rate' = FP/(TP+FP))
}
r1=18;n1=18;r2=6;n2=18
FPRPP.F(r1,r2,n1,n2,0.5)
p.value false positive rate
2.966259e-05 8.016273e-05
458:卵の名無しさん
18/04/11 21:24:11.18 gNeb3efk.net
# 非心カイ二乗分布
f9.22 <- function(k){
mui=runif(k)
X2.dash=0
f <- function(){
for(i in 1:k){
X2.dash=X2.dash+rnorm(1,mui[i])^2
}
return(X2.dash)
}
chisq.dash=replicate(10^3,f())
hist(chisq.dash,freq=FALSE,col=sample(colours(),1),xlab=quote(chi),main='')
ncp=sum(mui^2)
curve(dchisq(x,k,ncp),add=TRUE,lwd=2)
}
graphics.off()
x11() ; par(mfrow=c(3,3)) ; for(i in 1:9) f9.22(sample(3:10,1))
459:卵の名無しさん
18/04/12 00:50:35.65 GLkdI2IH.net
# S/σ2 ～ χ'(n-1,λ)
fp.142c <- function(n){
mui=rnorm(n)
sigma=runif(1)
mu=mean(mui)
xi=numeric(n)
f <- function(){
for(i in 1:n) xi[i]=rnorm(1,mui[i],sigma)
S=var(xi)*(n-1)
S/sigma^2
}
chisq.dash=replicate(10^3,f())
hist(chisq.dash,freq=FALSE,col=sample(colours(),1),xlab=quote(chi),main='')
ncp=sum((mui-mu)^2/sigma^2)
curve(dchisq(x,n,ncp),add=TRUE,lwd=2)
}
graphics.off()
fp.142c(7)
x11() ; par(mfrow=c(3,3)) ; for(i in 1:9) fp.142c(sample(3:10,1))
graphics.off()
460:卵の名無しさん
18/04/12 10:37:13.22 ZOlQ3FWn.net
URLﾘﾝｸ(imagizer.imageshack.com)
461:卵の名無しさん
18/04/12 18:48:37.23 GLkdI2IH.net
clt <- function(FUN,n=1000,k=10^4){
graphics.off()
x11()
par(mfrow=c(2,1))
hist(FUN(n),col='skyblue',freq=FALSE)
dat=replicate(k,FUN(n))
hist(apply(dat,1,sum),col=sample(colours(),1),freq=FALSE)
}
clt(FUN=function(x) rbeta(x,0.5,0.5))
clt(FUN=function(x) rexp(x))
clt(FUN=function(x) rgamma(x,1,1))
462:卵の名無しさん
18/04/15 01:25:53.84 2+v6Urh8.net
White=c(2.0477, 247, 1.8889)
Black=c(6.3772, 29, 1.4262)
Other=c(12.0832, 37, 3.7964)
# Totalc(3.6351 313 3.9770)
dat=rbind(White,Black,Other)
colnames(dat)=c("Mean","N","SD")
mean.G=sum(dat[,"Mean"]*dat[,"N"])/sum(dat[,"N"])
SS.bit=sum((dat[,"Mean"]-mean.G)^2*dat[,"N"])
SS.wit=sum(dat[,"SD"]^2*(dat[,"N"]-1))
df.bit=nrow(dat)-1
df.wit=sum(dat[,"N"]-1)
MS.bit=SS.bit/df.bit
MS.wit=SS.wit/df.wit
F.ratio=MS.bit/MS.wit
pf(F.ratio,df.bit,df.wit,lower.tail=FALSE)
(η2=(SS.bit)/(SS.bit+SS.wit))
Mean=dat[,'Mean'] ; N=dat[,'N'] ; SD=dat[,'SD']
sqrt(sum((N-1)*SD^2 + N*(Mean-sum(Mean*N)/sum(N))^2)/(sum(N)-1))
463:卵の名無しさん
18/04/15 10:56:43.94 2+v6Urh8.net
total.Mean <- function(Mean,N) sum(N*Mean)/sum(N)
total.SD <- function(Mean,N,SD) sqrt(sum((N-1)*SD^2 + N*(Mean-sum(Mean*N)/sum(N))^2)/(sum(N)-1))
(uniroot(function(x) total.Mean(c(60,65,x),c(20,5,95)) - 45,c(0,100))$root)
total.Mean(c(60,65,40.79),c(20,5,95))
(uniroot(function(x) total.SD(c(60,65,40.79),c(20,5,95),c(3,2,x)) - 10 , c(0,100))$root)
total.SD(c(60,65,40.79),c(20,5,95),c(3,2,6.13))
辞退者=c(60,20,3)
特待生=c(65, 5, 2)
裏口バカ=c(40.79, 95, 6.13)
# Totalc(45 120 10)
dat=rbind(辞退者,特待生,裏口バカ)
colnames(dat)=c("Mean","N","SD")
mean.G=sum(dat[,"Mean"]*dat[,"N"])/sum(dat[,"N"])
SS.bit=sum((dat[,"Mean"]-mean.G)^2*dat[,"N"])
SS.wit=sum(dat[,"SD"]^2*(dat[,"N"]-1))
df.bit=nrow(dat)-1
df.wit=sum(dat[,"N"]-1)
MS.bit=SS.bit/df.bit
MS.wit=SS.wit/df.wit
(F.ratio=MS.bit/MS.wit)
pf(F.ratio,df.bit,df.wit,lower.tail=FALSE)
(η2=(SS.bit)/(SS.bit+SS.wit)) # eta^2
> pf(F.ratio,df.bit,df.wit,lower.tail=FALSE)
[1] 2.789672e-30
> (η2=(SS.bit)/(SS.bit+SS.wit)) # eta^2
[1] 0.687539
>
464:卵の名無しさん
18/04/15 11:13:24.82 2+v6Urh8.net
辞退=scale(rnorm(20))*3+60
特待=scale(rnorm(5))*2+65
裏口=scale(rnorm(95))*6.1+40.79
score=c(辞退,特待,裏口)
group=c(rep('辞退',20),rep('特待',5),rep('裏口',95))
boxplot(score~group)
data=data.frame(score,group)
res=aov(score~group)
re=summary(res) ; re
str(re)
re2=re[[1]] ; re2
F.ratio=(re2$`Sum Sq`[1]/re2$`Df`[1])/(re2$`Sum Sq`[2]/re2$`Df`[2]) ; F.ratio
pf(F.ratio,re2$Df[1],re2$Df[2],lower.tail = FALSE)
465:卵の名無しさん
18/04/15 11:58:07.84 2+v6Urh8.net
##
Jonckheere<-function(L){
f<-function(A,B){
a<-length(A)
b<-length(B)
det<-0
for(i in 1:a){
for(j in 1:b){
det<- det+ifelse(A[i]==B[j],0.5,A[i]>B[j])
}
}
return(det)
}
g<-function(L){ #L=list(A1,,,,Aa),A1 > A2 > A3,..,> Aa : vector
a<-length(L)
comb<-combn(1:a,2)
con<-ncol(comb)
J=0
for(i in 1:con){
J<-J+f(L[[comb[1,i]]],L[[comb[2,i]]])
}
return(J)
}
466:卵の名無しさん
18/04/15 11:58:18.63 2+v6Urh8.net
J<-g(L)
a<-length(L)
n=double(a)
for(i in 1:a){
n[i]<-length(L[[i]])
}
N<-sum(n)
EJ <- (N^2-sum(n^2))/4
VJ <- (N^2*(2*N+3)-sum(n^2*(2*n+3)))/72
Z <- abs(J-EJ)/sqrt(VJ)
p.value<-pnorm(Z,lower.tail=FALSE)
return(p.value)
}
467:卵の名無しさん
18/04/15 16:37:48.33 2+v6Urh8.net
特待=scale(rnorm(5))*2+65
辞退=scale(rnorm(20))*3+60
裏口=scale(rnorm(95))*6.1+40.79
score=c(特待,辞退,裏口)
group=as.factor(c(rep(1,5),rep(2,20),rep(3,95)))
levels(group) <- c('特待','辞退','裏口')
boxplot(score~group)
data=data.frame(score,group)
res=aov(score~group)
re=summary(res) ; re
str(re)
re2=re[[1]] ; re2
F.ratio=(re2$`Sum Sq`[1]/re2$`Df`[1])/(re2$`Sum Sq`[2]/re2$`Df`[2]) ; F.ratio
pf(F.ratio,re2$Df[1],re2$Df[2],lower.tail = FALSE)
oneway.test(score~group)
kruskal.test(score~group,data=data)
Tk=TukeyHSD(aov(score~group,data=data)) ; Tk[[1]]
pairwise.t.test(score,group)
pairwise.t.test(score,group,p.adj='none')
pairwise.t.test(score,group,p.adj='bon')
pairwise.t.test(score,group,p.adj='bon',pool.sd = FALSE)
pairwise.wilcox.test(score,group)
#
library(PMCMR)
jonckheere.test(data$score,data$group,alternative = 'dec')
468:卵の名無しさん
18/04/16 08:27:01.84 9hc3WnzC.net
FPRP = Pr(H0|y)
＝ BF*PO/(BF*PO+1) BF = Pr(ｙ｜H0)/Pr(y | H1) : Bayes factor PO = π0/(1-π0) the prior odds of no association(true null hypothesis)
469:卵の名無しさん
18/04/16 11:44:56.69 j1Crkpch.net
αエラー、βエラーでのコストをCα,βCとすると
Pr(H0|y) < (Cβ/Cα)/(1+(Cβ/Cα))
A Bayesian Measure of the Probability of False Discovery
in Genetic Epidemiology Studies
470:卵の名無しさん
18/04/16 11:45:23.85 j1Crkpch.net
αエラー、βエラーでのコストをCα,Cβとすると
Pr(H0|y) < (Cβ/Cα)/(1+(Cβ/Cα))
が成立するとのこと。
A Bayesian Measure of the Probability of False Discovery
in Genetic Epidemiology Studies
471:卵の名無しさん
18/04/18 20:27:56.69 OaKjeXLn.net
f.nxP0<-function(x,n,P0,Ignore=FALSE){
P.hat=x/n
if((!Ignore)&(n*P.hat<5 | n*(1-P.hat)<5 | n*P0<5 | n*(1-P0)<5)){
warning('unfit for asymptotics')
return(NA)
}
u0= (P.hat-P0)/(sqrt(P0*(1-P0)/n))
p.value=pnorm(-abs(u0))
data.frame(u0,p.value)
}
> f.nxP0(x=3,n=100,P0=0.10)
[1] NA
Warning message:
In f.nxP0(x = 3, n = 100, P0 = 0.1) : unfit for asymptotics
> f.nxP0(3,100,0.10,Ignore=TRUE)
u0 p.value
1 -2.333333 0.009815329
472:卵の名無しさん
18/04/21 17:35:11.37 Hyq4QfSG.net
fisher.testと超幾何分布
x = 2 # of white balls drawn (<=m)
m = 6 # of white balls in urn
n = 4 # of black balls in urn
k = 5 # drawn
sum(dhyper(0:k,m,n,k)[dhyper(0:k,m,n,k)<=dhyper(x,m,n,k)])
fisher.test(matrix(c(x,k-x,m-x,n-k+x),ncol=2,byrow=TRUE))$p.value
473:卵の名無しさん
18/04/21 20:09:41.70 Hyq4QfSG.net
# URLﾘﾝｸ(oku.edu.mie-u.ac.jp)
?dhyper
x=2
wd=x # white drawn : x < td (= k)
wu=6 # white in urn : m
bu=4 #
474:blac in urn : n td=5 # total drawn : k wd2p <- function(wd) dhyper(wd,wu,bu,td) #white drawn to probability (probs=sapply(0:td,wd2p)) plot(0:td,probs,type='h') sum(probs[probs<=wd2p(x)]) fisher.test(matrix(c(wd,td-wd,wu-wd,bu-(td-wd)),ncol=2,byrow=TRUE))$p.value td=4 (probs=sapply(0:td,wd2p)) plot(0:td,probs,type='h') wd=1 sum(probs[probs<=wd2p(wd)]) fisher.test(matrix(c(wd,td-wd,wu-wd,bu-(td-wd)),ncol=2,byrow=TRUE))$p.value
475:卵の名無しさん
18/04/21 20:10:40.58 Hyq4QfSG.net
ex1 = matrix(c(6,12,12,5), nrow=2)
fisher.test(ex1)
exact2x2::fisher.exact(ex1)
476:卵の名無しさん
18/04/21 21:23:35.57 Hyq4QfSG.net
Event no Event
exposure a b
no exposure c d
OR=(a/c)/(b/d) = ad/bc ; ad=bc*OR
RR=a/(a+b) ÷ c/(c+d)=a(c+d)÷c(a+b)=(ac+ad)/(ac+bc)
=(ac+bc*OR)/(ac+bc)
=(b*OR+a)/(a+b)
=OR*b/(a+b)+a/(a+b)
=OR*( 1-a/(a+b) ) + a/(a+b)
=OR*(1-P1) + P1 # P1:暴露群でのイベント発生率=a/(a+b)
477:卵の名無しさん
18/04/21 21:25:20.08 Hyq4QfSG.net
Event no Event
exposure a b
no exposure c d
OR=(a/c)/(b/d) = ad/bc ; bc=ad/OR
RR=a/(a+b) ÷ c/(c+d)=a(c+d)÷c(a+b)=(ac+ad)/(ac+bc)
=(ac+ad)/(ac+ad/OR)
=(c+d) /(c + d/OR) # (c+d) で割る
=1 / {c/(c+d) + d/[(c+d)/OR]}
=OR/{OR*c/(c+d) + d/(c+d)} # ORを掛ける
=OR/{OR*c/(c+d) + (c+d)/(c+d)-c/(c+d)}
=OR/(OR*P0 + 1-P0) # P0:非暴露群でのイベント発生率=c/(c+d)
####
RR=OR*(1-P1) + P1
OR=(RR-P1)/(1-P1) # P1:暴露群でのイベント発生率=a/(a+b)
RR=OR/(OR*P0 + 1-P0)
OR=(1-P0)*RR/(1-P0*RR) # P0:非暴露群でのイベント発生率=c/(c+d)
####
OR2RRp0=function(OR,p0) OR/(OR*p0+1-p0) #p0:非暴露群でのイベント発生率
OR2RRp1=function(OR,p1) OR*(1-p1)+p1 #P1:暴露群でのイベント発生率
478:卵の名無しさん
18/04/21 21:32:57.14 Hyq4QfSG.net
##
RR=OR*(1-P1) + P1
OR=(RR-P1)/(1-P1) # P1:暴露群でのイベント発生率=a/(a+b)
P1=(OR-RR)/(OR-1)
RR=OR/(OR*P0 + 1-P0)
OR=(1-P0)*RR/(1-P0*RR) # P0:非暴露群でのイベント発生率=c/(c+d)
P0=(OR-RR)/(OR-1)/RR
RR2ORp1 <- function(RR,P1) (RR-P1)/(1-P1) # P1:暴露群でのイベント発生率=a/(a+b)
RR2ORp0 <- function(RR,P0) (1-P0)*RR/(1-P0*RR) # P0:非暴露群でのイベント発生率=c/(c+d)
R2p1 <- function(RR,OR) (OR-RR)/(OR-1)
R2p0 <- function(RR,OR) (OR-RR)/(OR-1)/RR
479:卵の名無しさん
18/04/22 08:01:09.15 TOqm84tJ.net
URLﾘﾝｸ(swada.net)
480:卵の名無しさん
18/04/22 23:30:16.86 TOqm84tJ.net
# 女子大生が全部で100人いるとする。
# 身長は170cmか，150cmのどちらかである。
# 身長170cmの人の胸囲は90cm, 身長150cmの人の胸囲は75cm．
# そして，胸囲90cmの人は確率0.9でクラス1へ，確率0.1でクラス2に振り分けられる.
# 一方で，胸囲75cmの人は確率0.8でクラス2へ，確率0.2でクラス1に振り分けられる.
# そして，胸囲が90cmと75cmの人はそれぞれ30,70人である。
# クラス1とクラス2の胸囲調整済み平均身長を求めよ。
# Class1
(170*30*0.9 + 150*70*0.2)/(30*0.9 + 70*0.2) # 27+14人
mean(c(rep(170,30*0.9),rep(150,70*0.2)))
# Class2
(170*30*0.1 + 150*70*0.8)/(30*0.1 + 70*0.8) # 3+56人
mean(c(rep(170,30*0.1),rep(150,70*0.8)))
# average
mean(c(rep(170,30),rep(150,70)))
# Inverse Probability Weighting(IPW)
# Class 1 adjusted
(170*1/0.9*27 + 150*1/0.2*14)/sum(1/0.9*27+1/0.2*14)
weighted.mean(c(rep(170,30*0.9),rep(150,70*0.2)),c(rep(1/0.9,30*0.9),rep(1/0.2,70*0.2)))
# Class 2 adjusted
(170*1/0.1*3 + 150*1/0.8*56)/sum(1/0.1*3+1/0.8*56)
weighted.mean(c(rep(170,30*0.1),rep(150,70*0.8)),c(rep(1/0.1,30*0.1),rep(1/0.8,70*0.8)))
481:卵の名無しさん
18/04/23 19:00:36.74 4SvCBV2R.net
Inverse Propensity score Weighting
IPW8 <- function(Y,Tr,ps) weighted.mean(Y,Tr/ps) # Y:dead or alive, Tr:treated or control, ps:propensity score
IPW8(Y,Tr,ps) - IPW8(Y,1-Tr,1-ps)
482:卵の名無しさん
18/04/24 12:13:34.39 6fgOc8MP.net
回帰分析による推定では，傾向スコアと，目的変数が線形な関係になる必要があるが，傾向スコア自体は（ロジスティック回帰による処置群に含まれる確率であるため）0-1の間の値をとるので，線形性を仮定するのは論理的におかしい
483:卵の名無しさん
18/04/24 22:47:05.61
484:5RekAiwJ.net
485:卵の名無しさん
18/04/25 10:26:10.03 5SphbBIi.net
# Lindner Center data on 996 PCI patients analyzed by Kereiakes et al. (2000)
library(MatchLinReg)
data("lindner")
head(lindner)
formula.stent=stent ~ abcix + height + female + diabetic + acutemi + ejecfrac + ves1proc-1
re.glm=glm(stent ~ abcix + height + female + diabetic + acutemi +ejecfrac + ves1proc-1,family=binomial,data=lindner)
summary(re.glm)
Epi::ROC(form=stent~.-1,data=lindner)$AUC
rms::lrm( stent ~ abcix + height + female + diabetic + acutemi +
ejecfrac + ves1proc, data=lindner)$stat['C']
ps=re.glm$fitted.values
Y=lindner$lifepres
Tr=lindner$stent
weighted.mean(Y,Tr/ps) - weighted.mean(Y,(1-Tr)/(1-ps))
486:卵の名無しさん
18/04/25 10:26:32.95 5SphbBIi.net
library(Matching)
m.out=Match(Y,Tr,ps)
summary(m.out)
## (Doubly Robust: DR)
n <- nrow(lindner)
y <- lindner['lifepres']
z <- lindner['stent']
lindner1 <- lindner[lindner['stent']==1,]
lindner0 <- lindner[lindner['stent']==0,]
model1=lm(formula=lifepres ~ abcix + height + female + diabetic
+ acutemi +ejecfrac + ves1proc, data=lindner1)
model0=lm(formula=lifepres ~ abcix + height + female + diabetic
+ acutemi +ejecfrac + ves1proc, data=lindner0)
fitted1 <- predict(model1, lindner)
fitted0 <- predict(model0, lindner)
dre1 <- (1/n)*sum(y+((z-ps)/ps)*(y-fitted1))
dre0 <- (1/n)*sum(((1-z)*y)/(1-ps)+(1-(1-z)/(1-ps))*fitted0)
c(treated=dre1,control=dre0,ATE=dre1-dre0)
#
487:卵の名無しさん
18/04/25 10:26:42.47 5SphbBIi.net
match0 <- function(Y,Tr,X,caliper=0.05){
interval = seq(0,1,by=caliper)
n=length(interval)
res = numeric(n)
for(i in 1:(n-1)){
temp0 = Y[ (Tr==0) & (interval[i]<X) & (X<interval[i+1]) ]
temp1 = Y[ (Tr==1) & (interval[i]<X) & (X<interval[i+1]) ]
res[i] = mean(temp1) - mean(temp0)
}
mean(res,na.rm=TRUE)
}
match0(Y=lindner$lifepres,Tr=lindner$stent,X=ps)
summary(Match(Y=lindner$lifepres,Tr=lindner$stent,X=ps, caliper=0.05))
MatchBalance(formula.stent,data=lindner)
488:卵の名無しさん
18/04/25 10:27:40.82 5SphbBIi.net
>>455
臨床統計のRスクリプトを公開するスレ。
489:卵の名無しさん
18/04/26 18:44:43.66 UTg/BRA4.net
ではHPVワクチン由来と断定できる重篤例に的を絞った質の高いエビデンスは得られるのでしょうか。
上記の相対リスク＝1.82を論文記載の対照群のリスク（36人/10万人）等を基に有意なエビデンスとして検出するには約19万人を対象とする臨床試験が必要です。
WHOにより検証された世界の質の高い臨床試験の総人数が7万3697人）であることを考えると，現時点で因果関係を高い質で問うのは困難であるとわかります
URLﾘﾝｸ(www.igaku-shoin.co.jp)
とあるが、
俺の計算では（Rのpackage、Epiで計算）
各群約46000人で
p=0.05
Relative Risk: 1.8200 1 3.3125
> Epi::twoby2(mat)
2 by 2 table analysis:
------------------------------------------------------
Outcome : Col 1
Comparing : Row 1 vs. Row 2
Col 1 Col 2 P(Col 1) 95% conf. interval
Row 1 30.1927 46051.50 7e-04 5e-04 9e-04
Row 2 16.5894 46065.11 4e-04 2e-04 6e-04
95% conf. interval
Relative Risk: 1.8200 1 3.3125
Sample Odds Ratio: 1.8205 1 3.3144
Probability difference: 0.0003 0 0.0006
Asymptotic P-value: 0.05
------------------------------------------------------
>
490:卵の名無しさん
18/04/26 18:58:06.22 UTg/BRA4.net
# URLﾘﾝｸ(www.igaku-shoin.co.jp)
# 相対リスク＝1.82を論文記載の対照群のリスク（36人/10万人）
# 1.82（95％信頼区間0.79-4.20）
rr=1.82
r0=36/10^5
r1=rr*r0
nn=seq(10^4,10^5,by=100)
rr2p.e <- function(n){
Epi::twoby2(matrix(c(n*r1,n*r0,n-n*r1,n-n*r0),ncol=2),print=FALSE)$p.value
}
plot(nn,sapply(nn,rr2p.e),type='l') ; abline(h=0.05, lty=3)
n=uniroot(function(n) rr2p.e(n)-0.05, c(10^4,10^5))$root
mat=matrix(c(n*r1,n*r0,n-n*r1,n-n*r0),ncol=2)
Epi::twoby2(mat)
491:卵の名無しさん
18/04/28 09:44:39.14 w/IREkFx.net
# 2 6 12 54 56+ 68 89 96 96 125+ 128+ 131+ 140+ 141+ 143 145+ 146 148+ 162+ 168 173+ 181+
time=c(2,6,12,54,56,68,89,96,96,125,128,131,140,141,143,145,146,148,162,168,173,181)
y=c(21/22,20/21,19/20,18/19,1,16/17,15/16,14/15,13/14,1,1,1,1,1,7/8,1,5/6,1,1,2/3,1,1)
survival=cumprod(y)
plot(time,survival,type='S',ylim=c(0,1))
censored=c(0,0,0,0,1,0,0,0,0,1,1,1,1,1,0,1,0,1,1,0,1,1)
event=1-censored
library(survival)
Surv(time,event)
plot(survfit(Surv(time,event)~1))
492:卵の名無しさん
18/04/28 11:58:57.59 w/IREkFx.net
カプランマイヤーのproduct limit法を理解していたらパッケージに頼らずにグラフも書けるな。
打ち切りポイントも描画できるように変更。
time=c(2,6,12,54,56,68,89,96,96,125,128,131,140,141,143,145,146,148,162,168,173,181)
y=c(21/22,20/21,19/20,18/19,1,16/17,15/16,14/15,13/14,1,1,1,1,1,7/8,1,5/6,1,1,2/3,1,1)
survival=cumprod(y)
censored=c(0,0,0,0,1,0,0,0,0,1,1,1,1,1,0,1,0,1,1,0,1,1)
plot(time,survival,type='s',ylim=c(0,1))
for(i in 1:22) if(censored[i]) points(time[i],survival [i],pch='+')
493:卵の名無しさん
18/04/30 11:29:32.46 T13xGV6f.net
# rule of thumb
# サンプルサイズと分散が同等な正規分布する集団からの無作為抽出で標本平均±標準誤差区間が重ならないときは母平均の有意差有無の判断はできない。
N=1000
f <- function(mx){
set.seed(123)
x1=rnorm(N,0 ,1)
x2=rnorm(N,mx,1)
t.test(x1,x2,var.equal = TRUE)$p.value
}
xx=seq(0,0.5,length=1000)
plot(xx,sapply(xx,f))
abline(h=0.05,lty=3)
uniroot(function(x) f(x)-0.05,c(0,1))$root
mx=0.0615 # 0.615以上で有意
g <- function(mx){
set.seed(123)
x1=rnorm(N,0 ,1)
x2=rnorm(N,mx,1)
n1=length(x1) ; n2=length(x2)
m1=mean(x1) ; m2=mean(x2)
s1=sd(x1) ; s2=sd(x2)
se1=s1/sqrt(n1) ; se2=s2/sqrt(n1)
(m1+se1) > (m2-se2) # overlap?
}
xx=seq(0,0.5,length=1000)
plot(xx,sapply(xx,g))
abline(h=0,lty=3)
uniroot(g,c(0,1))$root # non-overlap には0.0369以上
f(0.07) # 有意差あり
g(0.07) # overlapなし
f(0.05) # 有意差なし
g(0.05) # overlapなし
494:卵の名無しさん
18/04/30 12:22:04.07 T13xGV6f.net
# rule of thumb
# サンプルサイズと分散が同等な正規分布する集団からの無作為抽出で
# 標本平均±標準誤差区間が重なるときは母平均に有意差はない。
##
495:卵の名無しさん
18/04/30 14:24:19.80 T13xGV6f.net
# rule of thumb
# サンプルサイズと分散が同等な正規分布する集団からの無作為抽出で
# 標本平均±標準誤差区間が重なるときは母平均に有意差はない(p>0.05)。
#
x1=rnorm(n,m1,s)
x2=rnorm(n,m2,s) # m1 < m2
t.stat=(m1-m2)/(sqrt(2/n)*s) # < 0
=(m1-m2)*sqrt(n)/(sqrt(2)*s)
m1 + se > m2 -se # overlap
m1 - m2 > -2*se = -2*s/sqrt(n)
m1 - m2 = t.stat*(sqrt(2)*s)/sqrt(n) > -2*s/sqrt(n)
t.stat > -2/sqrt(2)=-sqrt(2)
df=198 # 例
curve(pt(x,df),-5,0) # x < 0で増加関数
pt(t.stat,df) > pt(-sqrt(2),df)
pt(-sqrt(2),df)
# dfを変化させる
curve(pt(-sqrt(2),x),200,ylim=c(0,0.5),xlab='df',ylab='p.value',lwd=2)
abline(h=0.05,lty=3)
pt(-sqrt(2),Inf) ; pnorm(-sqrt(2))
text(0,0.05,'0.05')
496:卵の名無しさん
18/04/30 14:24:33.69 T13xGV6f.net
# 平均値の95%信頼区間が重なった場合にはそれで有意差判定はできない。
set.seed(123)
x1=rnorm(100,0 ,1)
x2=rnorm(100,0.3 , 1)
t.test(x1,x2,var.equal = TRUE)$p.value # > 0.05
t.test(x1)$conf[2] > t.test(x2)$conf[1] # overlap TRUE
plot(NULL,NULL,ylim=c(0,0.05),xlim=c(-0.75,0.75),yaxt='n',ann=FALSE,bty='l')
segments(t.test(x1)$conf[1],0,t.test(x1)$conf[2],0,lwd=5,col=1)
segments(t.test(x2)$conf[1],0.01,t.test(x2)$conf[2],0.01,lwd=5,col=2)
set.seed(123)
x1=rnorm(100,0 ,1)
x2=rnorm(100,0.5 ,1)
t.test(x1,x2,var.equal = TRUE)$p.value # < 0.05
t.test(x1)$conf[2] > t.test(x2)$conf[1] # overlap TRUE
plot(NULL,NULL,ylim=c(0,0.05),xlim=c(-0.75,0.75),yaxt='n',ann=FALSE,bty='l')
segments(t.test(x1)$conf[1],0,t.test(x1)$conf[2],0,lwd=5,col=1)
segments(t.test(x2)$conf[1],0.01,t.test(x2)$conf[2],0.01,lwd=5,col=2)
497:卵の名無しさん
18/05/01 08:02:38.81 v/zhryAm.net
T.test2=function(n,dm,var1,var2){
SE12=sqrt((1/n+1/n)*((n-1)*var1+(n-1)*var2)/((n-1)+(n-1)))
T=dm/SE12
2*pt(abs(T),n-1+n-1,lower.tail = FALSE)
}
m1-m2=dm
n1=n2=n
T.test3=function(n,dm,var1,var2){
SE12=sqrt((2/n)*(var1+var2)/2))
T=dm/SE12
2*pt(abs(T),n-1+n-1,lower.tail = FALSE)
}
T.test4=function(x,n=1000,dm=1,var1=1){
SE=sqrt((2/n)*(var1+x)/2)
T=dm/SE
2*pt(abs(T),n-1+n-1,lower.tail = FALSE)
}
curve(T.test4(x),0,100,xlab='variance=x against variance=1',ylab='p.value')
498:卵の名無しさん
18/05/01 08:23:44.83 v/zhryAm.net
# t検定(生データなし,等分散不問)
Welch.test=function(n1,n2,m1,m2,var1,var2){
T=(m1-m2)/sqrt(var1/n1+var2/n2)
df=(var1/n1+var2/n2)^2 / (var1^2/n1^2/(n1-1)+var2^2/n2^2/(n2-1))
p.value=2*pt(abs(T),df,lower.tail = FALSE)
return(p.value)
}
n1=n2=n
m1-m2=dm
var2=x
Welch.test2=function(x,n=1000,dm=1,var1=1){
T=dm/sqrt(var1/n+x/n)
df=(var1/n+x/n)^2 / (var1^2/n^2/(n-1)+x^2/n^2/(n-1))
p.value=2*pt(abs(T),df,lower.tail = FALSE)
return(p.value)
}
curve(Welch.test2(x),0,100,xlab='variance=x against variance=1',ylab='p.value')
499:卵の名無しさん
18/05/01 09:37:21.50 UP3hBRO4.net
N=1000
n1=40
n2=60
SDR=10
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (10^4,f2.b())
hist (p.t)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (10^4,f2.bW())
hist (p.W)
500:卵の名無しさん
18/05/01 10:15:19.84 UP3hBRO4.net
N=1000
n1=40
n2=60
SDR=1
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (10^4,f2.b())
hist (p.t)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (10^4,f2.bW())
hist (p.W)
501:卵の名無しさん
18/05/01 10:21:41.36 UP3hBRO4.net
N=1000
n1=40
n2=60
SDR=1
dm=0.1
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR+dm
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (10^4,f2.b())
hist (p.t)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (10^4,f2.bW())
hist (p.W)
502:卵の名無しさん
18/05/01 12:51:04.77 UP3hBRO4.net
T.test4=function(var2,n1,n2,dm=0.5,var1=1){
SE12=sqrt((1/n1+1/n2)*((n1-1)*var1+(n2-1)*var2)/((n1-1)+(n2-1)))
T=dm/SE12
2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}
curve(T.test4(x,40,160,0.1),0,100,xlab='variance=x against variance=1',ylab='p.value')
503:卵の名無しさん
18/05/01 14:51:52.52 W1Ln0POm.net
# URLﾘﾝｸ(www.rips-irsp.com)
# Why Psychologists Should by Default Use Welch’s t-test
# Instead of Student’s t-test
par(mfrow=c(2,1))
N=1000
n1=10
n2=90
SDR=10
dm=1
k=10^4
set.seed(1234)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR+dm
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (k,f2.b())
hist (p.t,main='Student\'s t.test',col=sample(colours(),1))
mean(p.t < 0.05) # power(when dm!=0) or Type I error(when dm=0)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (k,f2.bW())
hist (p.W,main='Welch\'s t.test',col=sample(colours(),1))
mean(p.W < 0.05) # power(when dm!=0) or Type I error(when dm=0)
504:卵の名無しさん
18/05/01 14:57:05.14 W1Ln0POm.net
Why Psychologists Should by Default Use Welch’s t-test Instead of Student’s t-test
URLﾘﾝｸ(www.rips-irsp.com)
に触発されてシミュレーションの条件（サンプルサイズや分散比）を変えて追試してみた。
結果にびっくり
URLﾘﾝｸ(i.imgur.com)
あたかもド底辺シリツ医大卒を最高学府を履修したと呼ぶほどの差異に等しい。
Welchのrobustnessに再度感銘した。
505:卵の名無しさん
18/05/01 19:34:12.07 W1Ln0POm.net
par(mfrow=c(2,1))
N=1000
n1=50
n2=25
SDR=5
dm=0 # null hypothesis is true when dm==0
k=10^4
set.seed(123)
A=scale(rnorm(N))
B=scale(rnorm(N))*SDR+dm
f2.b <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=TRUE)$p.value
}
p.t=replicate (k,f2.b())
hist (p.t,main='Student\'s t.test',col=sample(colours(),1),
xlab=paste('n1 = ',n1,', sd1 = 1 ; n2 = ',n2,' , sd2 = ',SDR))
mean(p.t < 0.05) # power(when dm!=0) or Type I error(when dm=0)
f2.bW <- function(){
a=sample(A,n1)
b=sample(B,n2)
t.test(a,b, var=FALSE)$p.value
}
p.W=replicate (k,f2.bW())
hist (p.W,main='Welch\'s t.test',col=sample(colours(),1),
xlab=paste('n1 = ',n1,', sd1 = 1 ; n2 = ',n2,' , sd2 = ',SDR))
mean(p.W < 0.05) # power(when dm!=0) or Type I error(when dm=0)
506:卵の名無しさん
18/05/04 16:46:09.21 6BNcb2/S.net
x=log.nor=seq(-8.0,-3.0,by=0.5) ; n=length(x)
y=relax=c(2.6,10.5,15.8,21.1,36.8,57.9,73.7,89.5,94.7,100.0,100.0)
(re.nls=nls(y ~ Top/(1+10^((LogEC50-x)*HillSlope)),
start=c(Top=100,LogEC50=-5,HillSlope=1),algo='port'))
nls2R2(re.nls,y)
plot(x,y,bty='l',pch=19)
lines(x,predict(re.nls))
#
R2 <- function(dat){
n = nrow(dat)
sm = summary(
nls(dat$y ~ Top/(1+10^((LogEC50-dat$x)*HillSlope)),
start=c(Top=100,LogEC50=-5,HillSlope=0.5),algorithm = 'port')
)
SSalt = sum(sm$residuals^2)
SSnull = var(dat$y)*(n-1)
Rsq=1-SSalt/SSnull
Rsq
# k=sm$df[1]
# AdjRsq=1-(1-Rsq)*(n-1)/(n-k-1)
# data.frame(Rsq,AdjRsq)
}
xy
507:=data.frame(x,y) R2(xy) library(boot) re=boot(xy,function(data,indices) R2(data[indices,]),R=1000) re$t0 BEST::plotPost(re$t) ; HDInterval::hdi(re$t) quantile(re$t,c(.025,.975)) nls2R2(re.nls,y)
508:卵の名無しさん
18/05/04 17:55:57.07 6BNcb2/S.net
anova(lm)$Pr[1])
min( pf(summary(lm)$fstatistic[1],
summary(lm)$fstatistic[2],summary(lm)$fstatistic[2],lower=FALSE),
pf(summary(lm)$fstatistic[1],
summary(lm)$fstatistic[2],summary(lm)$fstatistic[2] )
509:卵の名無しさん
18/05/04 18:01:55.70 6BNcb2/S.net
lm = lm(formula,data)
lm2p <- function(lm){ # p_value=anova(lm)$Pr[1]
f=summary(lm)$fstatistic
min(
pf(f[1],f[2l,f[3],lower=FALSE),
pf(f[1],f[2l,f[3])
)
}
510:卵の名無しさん
18/05/04 20:28:04.76 6BNcb2/S.net
＃＃ ｈｔｔｐ：／／ｗｗｗ．ｐｕｂｌｉｃ．ａｓｕ．ｅｄｕ／～ｇａｓｗｅｅｔｅ／ｃｒｊ６０４／ｒｅａｄｉｎｇｓ／１９８３－Ｆｒｅｅｄｍａｎ％２０（Ｓｃｒｅｅｎｉｎｇ％２０Ｒｅｇｒｅｓｓｉｏｎ％２０Ｅｑｕａｔｉｏｎｓ）．ｐｄｆ
ｎｏｉｓｅ２ｓｉｇ ＜－ ｆｕｎｃｔｉｏｎ（ｓｅｅｄ，ｓｅｑｕｅｎｔｉａｌ＝ＦＡＬＳＥ）｛
ｓｅｔ．ｓｅｅｄ（ｓｅｅｄ）
Ｎｏｉｓｅ＝ｍａｔｒｉｘ（ｒｎｏｒｍ（１００＊５１），ｎｃｏｌ＝５１，ｎｒｏｗ＝１００）
＃ ｆｉｒｓｔ ｐａｓｓ
ｌｍ１＝ｌｍ（Ｎｏｉｓｅ［，５１］ ～ Ｎｏｉｓｅ［，－５１］） ＃ Ｎｏｉｓｅ［，５１］ ： 目的変数Ｙ
ｃａｔ（＇\ｎＲ２ （１ｓｔ） ＝ ＇，ｓｕｍｍａｒｙ（ｌｍ１）＄ｒ．ｓｑｕａｒｅｄ）
ｃａｔ（＇\ｎｐ．ｖａｌｕｅ （１ｓｔ） ＝ ＇，ａｎｏｖａ（ｌｍ１）＄Ｐｒ［１］）
ｃｏｅｆｓ１＝ｓｕｍｍａｒｙ（ｌｍ１）＄ｃｏｅｆ［，＇Ｐｒ（＞｜ｔ｜）＇］
ｌｅｎｇｔｈ（ｃｏｅｆｓ１）
ｃｏｅｆｓ１［１］
ｃｏｅｆｆｓ１＝ｃｏｅｆｓ１［－１］
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．２５ ＝＇，ｓｕｍ（ｃｏｅｆｆｓ１＜０．２５））
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．０５ ＝＇，ｓｕｍ（ｃｏｅｆｆｓ１＜０．０５））
ｉｎｄｘ２５＝ｗｈｉｃｈ（ｃｏｅｆｆｓ１＜０．２５）
511:卵の名無しさん
18/05/04 20:28:19.24 6BNcb2/S.net
＃ ｓｅｃｏｎｄ ｐａｓｓ
ｌｍ２＝ｌｍ（Ｎｏｉｓｅ［，５１］ ～ Ｎｏｉｓｅ［，ｉｎｄｘ２５］）
ｃａｔ（＇\ｎ\ｎＲ２ （２ｎｄ） ＝ ＇，ｓｕｍｍａｒｙ（ｌｍ２）＄ｒ．ｓｑｕａｒｅｄ）
ｃａｔ（＇\ｎｐ．ｖａｌｕｅ（２ｎｄ） ＝ ＇，ａｎｏｖａ（ｌｍ２）＄Ｐｒ［１］）
ｃｏｅｆｓ２＝ｓｕｍｍａｒｙ（ｌｍ２）＄ｃｏｅｆ［，＇Ｐｒ（＞｜ｔ｜）＇］
ｌｅｎｇｔｈ（ｃｏｅｆｓ２）
ｃｏｅｆｓ２［１］
ｃｏｅｆｆｓ２＝ｃｏｅｆｓ２［－１］
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．２５ （２ｎｄ） ＝＇，ｓｕｍ（ｃｏｅｆｆｓ２＜０．２５））
ｃａｔ（＇\ｎｃｏｅｆ ＜ ０．０５ （２ｎｄ） ＝＇，ｓｕｍ（ｃｏｅｆｆｓ２＜０．０５））
ｃａｔ（＇\ｎ\ｎ＇）
ｉｆ（ｓｅｑｕｅｎｔｉａｌ）｛
ｃａｔ（＇Ｈｉｔ Ｒｅｔｕｒｎ Ｋｅｙ ｉｎ ｃｏｎｓｏｌｅ ｗｉｎｄｏｗ＇）
ｎｏ＿ｓａｖｅ ＜－ ｓｃａｎ（ｎ＝１， ｗｈａｔ＝ｃｈａｒａｃｔｅｒ（）， ｑｕｉｅｔ＝ＴＲＵＥ）
｝
｝
ｎｏｉｓｅ２ｓｉｇ（１）
ｆｏｒ（ｉ ｉｎ ２：５） ｎｏｉｓｅ２ｓｉｇ（ｉ，ｓｅｑ＝ＴＲＵＥ）
512:卵の名無しさん
18/05/07 08:41:59.73 /y3oEDZQ.net
UG.Purge <- function(alt='two.sided'){
pv=numeric(10)
for(j in 1:10){
L=list()
tmp=UG[-j,]
for(i in 1:6) L=append(L,list(tmp[,i]))
pv[j]=jonckheere(L,cat=FALSE,alternative = alt)
}
return(pv)
}
which(UG.Purge('two') < 0.05)
which(UG.Purge('increasing') < 0.05)
513:卵の名無しさん
18/05/08 17:13:30.89 4Ru4wKk7.net
＃ ｈｔｔｐ：／／ｆｉｌｅ．ｓｃｒａｔｃｈｈｉｔ．ｐａｚｒｕ．ｃｏｍ／ＴｋｙＳｃｈＢｎｆ．ｔｘｔ
＃＃ ｆｉｔｔｅｄ ｉｓ ａ ｓｅｑｕｅｎｃｅ ｏｆ ｍｅａｎｓ
＃＃ ｎｉｓ ｉｓ ａ ｃｏｒｒｅｓｐｏｎｄｉｎｇ ｓｅｑｕｅｎｃｅ ｏｆ ｓａｍｐｌｅ ｓｉｚｅｓ ｆｏｒ ｅａｃｈ ｍｅａｎ
＃＃ ｄｆ ｉｓ ｔｈｅ ｒｅｓｉｄｕａｌ ｄｆ ｆｒｏｍ ｔｈｅ ＡＮＯＶＡ ｔａｂｌｅ
＃＃ ＭＳＥ ＝ ｍｅａｎ ｓｑｕａｒｅｄ ｅｒｒｏｒ ｆｒｏｍ ｔｈｅ ＡＮＯＶＡ ｔａｂｌｅ
＃＃ ｃｏｎｆ．ｌｅｖｅｌ ｉｓ ｔｈｅ ｆａｍｉｌｙ－ｗｉｓｅ ｃｏｎｆｉｄｅｎｃｅ ｌｅｖｅｌ， ｄｅｆａｕｌｔｓ ｔｏ ．９５
ｐａｉｒｗｉｓｅ．ｓｃｈｅｆｆｅ．ｔ�
514:�ｓｔ ＜－ ｆｕｎｃｔｉｏｎ（ｒｅ．ａｏｖ， ｃｏｎｆ．ｌｅｖｅｌ ＝ ０．９５）｛ ｍｏｄｅｌ ＝ ｒｅ．ａｏｖ＄ｍｏｄｅｌ ｃｏｌｎａｍｅｓ（ｍｏｄｅｌ） ＝ ｃ（＇ｓｃｏｒｅ＇，＇ｇｒｏｕｐ＇） ｆｉｔｔｅｄ ＝ ｗｉｔｈ（ｍｏｄｅｌ，ｔａｐｐｌｙ（ｓｃｏｒｅ，ｇｒｏｕｐ，ｍｅａｎ）） ｎｉｓ ＝ ｗｉｔｈ（ｍｏｄｅｌ，ｔａｐｐｌｙ（ｓｃｏｒｅ，ｇｒｏｕｐ，ｌｅｎｇｔｈ）） ｄｆ ＝ ｒｅ．ａｏｖ＄ｄｆ．ｒｅｓｉｄｕａｌ ＭＳＥ ＝ ｓｕｍｍａｒｙ（ｒｅ．ａｏｖ）［［１］］［＇Ｒｅｓｉｄｕａｌｓ＇，＇Ｍｅａｎ Ｓｑ＇］ ｒ ＝ ｌｅｎｇｔｈ（ｆｉｔｔｅｄ） ｐａｉｒｓ ＝ ｃｏｍｂｎ（１：ｒ，２） ｄｉｆｆｓ ＝ ｆｉｔｔｅｄ［ｐａｉｒｓ［１，］］ － ｆｉｔｔｅｄ［ｐａｉｒｓ［２，］］
515:卵の名無しさん
18/05/08 17:13:57.38 4Ru4wKk7.net
Ｔ ＝ ｓｑｒｔ（（ｒ－１）＊ｑｆ（ｃｏｎｆ．ｌｅｖｅｌ，ｒ－１，ｄｆ））
ｈｗｉｄｔｈｓ ＝ Ｔ＊ｓｑｒｔ（ＭＳＥ＊（１／ｎｉｓ［ｐａｉｒｓ［１，］］ ＋ １／ｎｉｓ［ｐａｉｒｓ［２，］］））
ｆｖｓ ＝ （ｄｉｆｆｓ＾２）／（ＭＳＥ＊（１／ｎｉｓ［ｐａｉｒｓ［１，］］ ＋ １／ｎｉｓ［ｐａｉｒｓ［２，］］））／（ｒ－１）
ｐｖｓ ＝ １－ｐｆ（ｆｖｓ， ｒ－１， ｄｆ）
ｖａｌ ＝ ｃｂｉｎｄ（ｄｉｆｆｓ － ｈｗｉｄｔｈｓ， ｄｉｆｆｓ， ｄｉｆｆｓ ＋ ｈｗｉｄｔｈｓ， ｆｖｓ， ｒ－１， ｄｆ， ｐｖｓ）
ｄｉｍｎａｍｅｓ（ｖａｌ） ＝ ｌｉｓｔ（ｐａｓｔｅ（＂ｓｕｂｊｅｃｔ＂，ｐａｉｒｓ［１，］，＂ － ｓｕｂｊｅｃｔ＂， ｐａｉｒｓ［２，］，
ｓｅｐ ＝ ＂＂）， ｃ（＂Ｌｏｗｅｒ＂， ＂Ｄｉｆｆ＂，＂Ｕｐｐｅｒ＂， ＂Ｆ Ｖａｌｕｅ＂， ＂ｎｍｄｆ＂， ＂ｄｎｄｆ＂， ＂Ｐｒ（＞Ｆ）＂））
Ｓｃ ＝ ａｓ．ｍａｔｒｉｘ（ｖａｌ）
ｐ．Ｓｃ ＝ Ｓｃ［，＇Ｐｒ（＞Ｆ）＇］
ｎ ＝ ｒｅ．ａｏｖ＄ｒａｎｋ － １
ｍａｔ．Ｓｃ ＝ ｍａｔｒｉｘ（ｒｅｐ（ＮＡ，ｎ＊ｎ），ｎ）
ｆｏｒ（ｋ ｉｎ １：ｎ） ｍａｔ．Ｓｃ［，ｋ］ ＝ ｃ（ｒｅｐ（ＮＡ，ｋ－１）， ｐ．Ｓｃ［（（ｋ－１）＊（ｎ＋１－ｋ／２）＋１）：（（ｋ－１）＊（ｎ＋１－ｋ／２）＋１＋ｎ－ｋ）］）
ｃｏｌｎａｍｅｓ（ｍａｔ．Ｓｃ） ＝ ｍｏｄｅｌ［［２］］［１：ｎ］
ｒｏｗｎａｍｅｓ（ｍａｔ．Ｓｃ） ＝ ｍｏｄｅｌ［［２］］［２：（ｎ＋１）］
ｄ ＝ ｒｏｕｎｄ（ｍａｔ．Ｓｃ，５）
ｄ［ｉｓ．ｎａ（ｄ）］ ＝ ＇－＇
ｐｒｉｎｔ（ｄ，ｑｕｏｔｅ ＝ ＦＡＬＳＥ）
ｉｎｖｉｓｉｂｌｅ（ｖａｌ）
｝
516:卵の名無しさん
18/05/09 14:26:52.13 L5gMx8M1.net
# disese, No disease
# positive test TP(A) FP(B)
# negative test FN(C) TN(D)
PV <- function(prevalence=c(10^-4,1),sn=0.80,sp=0.95,...){
plot(NULL,NULL,xlim=prevalence,ylim=c(0,1),xlab='Prevalence',
ylab='Predicative Value(Test Credibility)',type='n',bty='l',
main=paste0('感度 = ',sn, ' 特異度 = ',sp), ...)
legend('center',bty='n',lwd=2,lty=1:3,legend=c('陽性適中率','陰性適中率','偽陽性率'),
col=c('red','darkgreen','brown'))
ppv<-function(prevalence,sensitivity=sn,specificity=sp){
PPV=prevalence*sensitivity/
(prevalence*sensitivity+(1-prevalence)*(1-specificity))
return(PPV)
}
curve(ppv(x),lty=1,lwd=2,col='red', add=TRUE)
npv<-function(prevalence,sensitivity=sn,specificity=sp){
NPV=(1-prevalence)*specificity/
( (1-prevalence)*specificity + prevalence*(1-sensitivity) )
return(NPV)
}
curve(npv(x),lty=2,lwd=2,col='darkgreen',add=TRUE)
false_negative_rate <- function(prevalence,sensitivity=sn,specificity=sp){
FNR <- prevalence*(1-sensitivity)/
( (1-prevalence)*specificity + prevalence*(1-sensitivity) )
return(FNR)
}
curve(false_negative_rate(x), lty=3,lwd=2,col='brown',add=TRUE)
}
517:卵の名無しさん
18/05/11 20:34:40.05 MdgPLO2C.net
zodiac=c('Aquarius','Aries','Cancer','Capricorn','Gemini','Leo','Libra','Pisces','Sagittarius','Scorpio','Taurus','Virgo')
N=c(856301,888348,917553,844635,937615,903009,897503,893332,846813,850128,918512,921196)
X=c(1433,1476,1496,1343,1553,1497,1350,1522,1277,1297,1534,1445)
iP=which(zodiac=='Pisces')
p.Pisces=prop.test(c(X[iP],sum(X[-iP])),c(N[iP],sum(N[-iP])))$p.value
p0=sum(X)/sum(N)
se
518:t.seed(16062006) p.max=numeric() for(i in 1:10^4){ xi=rbinom(length(N),N,p0) I=which.max(xi/N) re=prop.test(c(xi[I],sum(xi[-I])),c(N[I],sum(N[-I]))) p.max[i]=re$p.value } mean(p.max <= p.Pisces) cat("P-value for Pisces and CHF: ", p.max, file="Pisces.txt",fill=TRUE, append=TRUE)
519:卵の名無しさん
18/05/13 14:56:25.30 lYhicRHD.net
TRAP9p <- function(seed=NULL){
set.seed(seed)
A0=rnorm(10,50,10)
A1=A0+rnorm(10,0,5)+5
B0=rnorm(10,45,10)
B1=B0+rnorm(10,0,5)+5
A0=round(A0)
B0=round(B0)
A1=round(A1)
B1=round(B1)
a=t.test(A0,A1,paired=TRUE)$p. # 外部有意差なし
b=t.test(B0,B1,paired=TRUE)$p. # 内部有意差あり
ab=t.test(A1-A0,B1-B0)$p. # 加点差有意差なし
a0b0=t.test(A0,B0)$p.# base差なし
a1b1=t.test(A1,B1)$p.# final差
list(p=c(外部=a,内部=b,差異=ab,base差=a0b0,final差=a1b1),A0=A0,B0=B0,A1=A1,B1=B1)
}
TRAP9p()$p
rep=replicate(1000,TRAP9p()$p)
mean(rep['外部',]>0.05 & rep['内部',]<0.05 & rep['差異',]>0.05
& re['base差',] > 0.05 & re['final差',]<0.05)
seeds=c(4799, 10638, 12173, 12908, 13671, 17145, 18955)
re=sapply(seeds,function(x) TRAP9p(x)$p) # wait several ten seconds
idx=which(re['外部',]>0.05 & re['内部',]<0.05 & re['差異',]>0.05
& re['base差',] > 0.05 & re['final差',]<0.05)
idx ; re[,idx]
Tx=as.data.frame(TRAP9p(4799)[2:5]) # 4799 10638 12173 12908 13671 17145 18955
U=with(Tx,data.frame(外部.無=A0,外部.有=A1,内部.無=B0,内部.有=B1))