臨床統計もおもしろいですよ、その3
内科認定医受験の最低限の知識、 製薬会社の示してくる臨床データ、 論文の考察、 論文を書くときの正当性、 というのが、臨床統計の今までの目的の大きい部分でしたが、 AI=機械学習の基本も、結局は統計学と確率に支配されます。 そういう雑多な話をするスレです。 ※前スレ 臨床統計もおもしろいですよ、その1 https://egg.2ch.net/test/read.cgi/hosp/1493809494/ 臨床統計もおもしろいですよ、その2 https://egg.5ch.net/test/read.cgi/hosp/1540905566/
中央値比較での統計悪用シミュレーションプログラム sim <- function(){ A=c(sort(sample(1:4,17,rep=T)),4,sort(sample(5:50,17,rep=T))) C=c(sort(sample(1:11,23,rep=T)),11,sort(sample(11:20,23,rep=T))) print(summary(A)) print(summary(C)) par(mfrow=c(2,1)) plot(table(A),bty='l',xlim=c(1,50),col='red') plot(table(C),bty='l',xlim=c(1,50),col='blue') } sim() >>64 ド底辺頭脳はサイコロを6回投げたら1の目が1回でているとう風な計算しかできんみたいだな。 信頼区間とか確率分布という概念での思考ができなんよね。 周辺度数固定して二項分布で乱数発生させてχ二乗検定したときのp値の分布は以下のようになる。 https://i.imgur.com/5djScD1.png https://i.imgur.com/RsFr7o4.png 危険率1%で有効といえる確率はほぼ5割であることがわかる。 >>60 エントリーに5以下の数値があるとχ二乗検定だとYatesの補正を使うけど、この程度の標本数は臨床では普通だよ。 中国のアビガンの予備試験も標本数が実薬群ととコントロール群を合わせて80だぞ。 足りないのはオマエの学力だよ。 >>70 まあ 学園カースト最下位のFランさんが必死ですわね 昨日は先輩方からも国試合格のご報告がたくさんありましたが Fランさんのお知り合いの方たちはいかがですかぁ?W by 都内Sラン女子高生 >>71 答えられない低学力を隠そうと必死だな。 どうやって解くかわからんの? >>72 わかりました 犯人はアナタですね Fランさん ようやく完成した # 球面上の一様分布 vertex <- function(r=1){ x=runif(1,-1,1) # x ~ 一様分布[-1,1] phi=runif(1,-pi,pi) # φ ~ 一様分布[-π,π] y=sqrt(1-x^2)*cos(phi) # √(1-x^2)*cos(φ) z=sqrt(1-x^2)*sin(phi) # √(1-x^2)*sin(φ) r*c(x,y,z) } https://i.imgur.com/xX0mTim.png https://i.imgur.com/H7hs9w8.png # 凸包四面体の体積 sim <- function(r=1,print=F){ v4=replicate(4,vertex()) # 4点の直交座標 if(print) print(v4) abs(det(v4[,2:4]-v4[,1]))/6 # 四面体の体積 } sim(print=T) k=1e5 tetra=replicate(k,sim()) # k回のシミュレーション mean(tetra) summary(tetra) BEST::plotPost(tetra) # 単位球体内部一様分布の4点で四面体をつくるとき rm(list=ls()) vertex <- function(){ x=runif(1,-1,1) # x ~ 一様分布[-1,1] phi=runif(1,-pi,pi) # φ ~ 一様分布[-π,π] r=runif(1,0,1) # r ~ 一様分布[0,1] x=r^(1/3)*x # r^(1/3)*x y=r^(1/3)*sqrt(1-x^2)*cos(phi) # r^(1/3)√(1-x^2)*cos(φ) z=r^(1/3)*sqrt(1-x^2)*sin(phi) # r^(1/3)√(1-x^2)*sin(φ) c(x,y,z) } vtx=replicate(5000,vertex()) par(mfrow=c(3,1)) x=vtx[1,] ; hist(x,col='pink') y=vtx[2,] ; hist(y,col='orange') z=vtx[3,] ; hist(z,col='darkgreen') rgl::plot3d(x,y,z, col="skyblue") par(mfrow=c(1,1)) # 四面体の体積 sim <- function(r=1,print=F){ v4=replicate(4,vertex()) # 4点の直交座標 if(print) print(v4) abs(det(v4[,2:4]-v4[,1]))/6 # 四面体の体積 } sim(print=T) k=1e5 tetra=replicate(k,sim()) # k回のシミュレーション mean(tetra) summary(tetra) BEST::plotPost(tetra) core i7 Memory 16G はシミュレーションが捗って( ・∀・)イイ!! 臨床医学は確率事象を扱う稼業なので統計処理は必須。 これができない医師はアホといえる。 某国の新型コロナ感染症の有病率を0.1、 PCR検査の感度を0.7 特異度を0.9とする 検査陽性陰性の人を無作為に100人ずつ集めて、カレーを頻回に食べているかを調査した結果が 以下の通りであった。 検査陽性 検査陰性 カレー頻食 a(=36) b(=60) カレー稀食 c(=64) d(=40) カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか? Parameter Total (n=83) Severe/critica l Group (n=25) Ordinary Group (n=58) P Value Number of lobes 5 (4-5) 5 (5-5) 5 (3-5) 0.003 median (interquartile range, IQR) library(coin) # coin::wilcox_test(c(Ea,La) ~ factor(c(rep('Ea',n1),rep('La',n2)))) f4 <- function(x){ # filer p.value=0.003 # pv=t.test(sc[x[1],],or[x[2],])$p.value SC=sc[x[1],] OR=or[x[2],] w=coin::wilcox_test(c(SC,OR) ~ factor(c(rep('SC',25),rep('OR',58))),distribution="exact") pv=coin::pvalue(w) if(round(pv,3)==0.003) return(c(x[1],x[2])) else(return(c(0,0))) } PCR <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){ # n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画 prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定 PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算 adjtd.prev = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算 # NPVでも同じ結果 # NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算 # adjtd.prev = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算 mean=mean(adjtd.prev) # 期待値 median=median(adjtd.prev) # 中央値 mode=density(adjtd.prev)$x[which.max(density(adjtd.prev)$y)] # 最頻値 LU=unlist(HDInterval::hdi(adjtd.prev))[1:2] # 95%信頼区間 re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値 if(print){ # ヒストグラム描画 par(mfrow=c(2,2)) hist(prev,freq=F,xlim=c(0,0.5), main='prevalence') ; lines(density(prev)) hist(PPV,freq=F,xlim=c(0,0.5),main='PPV') ; lines(density(PPV)) hist(adjtd.prev,freq=F,xlim=c(0,0.5),main='adjusted prevalence ') ; lines(density(adjtd.prev)) BEST::plotPost(adjtd.prev) ; lines(density(adjtd.prev),col='skyblue') par(mfrow=c(1,1)) } print(round(re,4)) # 概算値表示 invisible(list(re,adjtd.prev)) # 代表値と事後有病率を非表示で返す } PCR(100,1) PCR(1000,10) PCR(10000,100) PCR(100000,1000) # stanでMCMCさせることにした。 data{ // corona.stan int n; // 1000 int x; // 10 real<lower=0,upper=1> sen; // 0.7 real<lower=0,upper=1> spc; // 0.9 } parameters{ real<lower=0,upper=1> prev; // prevalence } transformed parameters{ real<lower=0,upper=1> p; p = prev*sen + (1-prev)*(1-spc) ; // probability of positive test result } model{ x ~ binomial(n,p); prev ~ beta(1,1); // prev ~ uniform(0,1) } stanはコンパイルに時間がかかって面倒。 jagsの方が早くて( ・∀・)イイ!! rm(list=ls()) graphics.off() options(scipen = 4) options(digits=5) library(rjags) N=1000 ; X=10 dataList=list(n=N,x=X,sen=0.7,spc=0.9) modelstring=paste0(' model { prev ~ dbeta(1,1) p <- prev*sen + (1-prev)*(1-spc) x ~ dbin(p,n) } ') writeLines(modelstring,'TEMPmodel.txt') jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList) update(jagsModel) codaSamples = coda.samples( jagsModel , variable=c("prev","p"), n.iter=50000 ) js=as.matrix(codaSamples) BEST::plotPost(js[,'prev']) round(HDInterval::hdi(js[,'prev'])[1:2],7) 何の責任もとれないド底辺ほど 他人の責任にうるさく自分の権利ばかりしつこく要求しがちですわね コンプ薬屋のことですけど Fランのアナタも同類ですわよね 早く呼んでくればぁ by 都内Sラン女子高生 >>1 よそのスレで騒いでるみたいですけど Fランのアナタには知性があるとは 医師の方々は考えてはいないようですよ by 都内Sラン女子高生 >>1 国試に行き詰まって統計遊びですか こういうときハゲワロスって言うんでしたっけ? by ドSなSラン女子高生 タバコを吸わないものは自室に禁煙の張り紙など貼らない訳ですが 裏口裏口と盛んにワメキ回るアナタには 身に覚えがあると言うことでよろしいですねW 容疑者はオマエだぁっ(コナ○口調で) by ドSなSラン女子高生 アナタの期待値は100発0中だぁっ ち○ぴょろすぽ○ん by ドSなSラン女子高生 自分で作った統計ごっこのスレがここにあるのに なーんでよそのスレに迷惑をかけに行くのかなん やっぱりFラン脳のゆえんなのかしらねんのねん by ドSなSラン女子高生 >>88 これやってみ、ド底辺頭脳のオマエには無理だろうな。 24時間待ってやるぞ。 ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか? 正確な感度も特異度もわからないから有病率は計算できないという主張が数学板でされたから、 感度・特異度を確率変数としてベイズ階層モデルを組めば計算できること投稿しておいた。 最頻値と標準偏差からベータ分布のパラメータを計算するのに連立方程式を解く方が面倒だった。 きゃー、統計ごっこのクルクルパー来たーっW (キモワロス) ド底辺仲間のコンプ薬屋はまだですか? by ドSなSラン女子高生 >>92 >>88 これやってみ、ド底辺頭脳のオマエには無理だろうな。 24時間待ってやるぞ。 ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか? >>93 国試に合格できなかった還暦ジジイって本当ですかぁ? アチコチで反日レスしてるみたいですねー by ドSなSラン女子高生 >>95 Fラン高の宿題はFランのアナタ自身でされてくださいね Sランのアタシはアナタと違ってヒマではございませんのでW by ドSなSラン女子高生 >>96 ドツボ問題もできないドアホ。 これやってみ、ド底辺頭脳のオマエには無理だろうな。 24時間待ってやるぞ。 ドツボ13は100発0中 ドツボ14は10発0中 ドツボ15は1発0中 とする。 各々10000発撃ったときドツボの命中数の期待値はいくらか? 95人調べて47人陽性。感度50-70% 特異度90%前後として この集団の有病率はstanでMCMCして計算すると mean lower upper 0.7558614 0.5214634 0.9999897 感度60% 特異度90%に固定すると > pn2ip(95,47) infected prevalence 75.0000000 0.7894737 麻薬取締法違反 札幌の薬剤師と法人を略式起訴 病院で管理する医療用麻薬の数量について道に虚偽の届け出をしたとして、 札幌区検は19日までに麻薬取締法違反の罪で、札幌ひばりが丘病院(札幌 市厚別区)に勤務していた30代の薬剤師の男と、同病院を運営する医療法人 潤和会(同区)を略式起訴した。札幌簡裁がそれぞれ20万円以下の罰金刑 を言い渡す見通し。略式命令請求によると、男は2015年11月、道に対し、 院内で使用した麻薬の使用量を偽って届け出るなどしたとされる。潤和会には 法人も罰する同法の両罰規定を適用した。 出典:北海道新聞 平成30年10月19日付 晋型コロナ肺炎に感度0.9,特異度0.9の迅速検査が開発されたと仮定する。 日本人1億2595万人からX人を無作為抽出して有病率を推定したい。 有病率の99%信頼区間幅を1%以内で検定したい。 何人を抽出すれば十分といえるか? ss2sbj <- function(SEN,SPC,range=0.01,CONF=0.99,print=FALSE){ pr2sbj <- function(p0,sen=SEN,spc=SPC,R=range,conf.level=CONF){ z = qnorm(1-(1-conf.level)/2) p = p0*sen+(1-p0)*spc n = 4*z^2*p*(1-p)/R^2 return(n) } if(print){ prevalence=seq(0,1,by=0.001) plot(prevalence,sapply(prevalence,pr2sbj), ylab='subjects',bty='l',type='l',lwd=2)} optimize(function(x) pr2sbj(x,SEN,SPC), c(0,1) ,maximum = TRUE) } "n(=100)人の中から無作為にm(=10)人選んだらその中に少なくとも一人の感染者がいた。 全体で何人の感染者がいるかの期待値と95%信頼区間を求めよ。" rm(list=ls()) fn <- function(n,m,print=FALSE){ library(gmp) # 離散量で計算 x=0:n # 感染者数:x, 非感染数:n-x pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (m人全員非感染の確率) pdf=pmf/sum(pmf) # 確率密度関数化して Ev=sum(x*pdf) # 期待値を計算 if(print) print(Ev) E=as.numeric(Ev) # 信頼区間計算計算のために連続量とする y=seq(0,n,by=0.05) # infected p=1- choose(n-y,m)/choose(n,m) # Pr[at least one infected among m] d=p/sum(p) # pmf -> pdf (e=sum(y*d)) # mean of probability if(print){plot(y,d,bty='l',xlab='infected',ylab='density',type='l') print(quantile(n*p,prob=c(0.025,0.05,0.5,0.95,0.975))) print(HDInterval::hdi(n*p)[1:2])} return(E) } fn(100,5,print=T) n=100 m=1:n plot(m,sapply(m,function(x) fn(n,x)), type='p',pch=19,bty='l',ylim=c(0,100), xlab='有名人',ylab='期待値',main='総人口:100人') " n(=10)人の中から無作為にm(=2)人選んだらその中に少なくとも一人の感染者がいた。 全体で何人の感染者がいるかの期待値を求めよ。 " f <- function(p){ # p:感染確率 p1=2*p*(1-p) # 一人だけ感染確率 p2=p^2 # 二人とも感染確率 (1*p1+2*p2)/(p1+p2) # 感染人数の期待値 } curve(f(x),bty='l',lwd=2,xlab='p') f(2/3) # 有名人が感染 " 長時間必要なのて実行非推奨,メモリ16G以下はクラッシュする library(gmp) m=18200 # 有名人の数(桜を見る会参加人数) n=1.268e5 # 日本の人口 x=0:n # 感染者数:x, 非感染数:n-x pmf=try(1- chooseZ(n-x,m)/chooseZ(n,m)) # 1 - (m人全員非感染の確率) pdf=try(pmf/sum(pmf)) # 確率密度関数化して (E=try(sum(x*pdf))) # 期待値を計算 as.numeric(E) # E=63463.27 (m=1000) , E=1154070201/18202=63403.48(m=1.268e5) " fn(100,20,print=T) data.frame(有名人=1:10,期待値=sapply(1:10,function(x) fn(100,x)$mean)) data.frame(有名人=1:10*10,期待値=sapply(1:10*10,function(x) fn(100,x)$mean)) fn4 <- function(n,m,k,quick=TRUE,print=FALSE){ # Pr[B|Ax] n人中x人感染者がいるときにm人中k人感染している確率 library(gmp) x=0:n if(quick){ Ev=sum(x*chooseZ(x,k))/sum(chooseZ(x,k)) }else{ pmf=chooseZ(x,k)/chooseZ(n,m) pdf=pmf/sum(pmf) Ev=sum(x*pdf)} # Σ(x*(xCk/nCm))/Σ(xCk/nCm) = Σ(x*(xCk))/Σ(xCk) if(print) print(Ev) E=as.numeric(Ev) return(E) } fn4=Vectorize(fn4) fn4(100,10,0:10,quick=F) 休校続きでお勉強くらいしかすることがございません 今思うのは反日野党が政権をとってたら 世界はモノスゴクひどいことになっているであろうこと by 都内Sラン女子高生 今にして思えば菅直人は実に優秀だったな。 事故当日に原子力緊急事態宣言を発令、命をかけて原発に突入して手動ベントを約束させる、 自衛隊10万に派遣決定、 全原発停止、 東電の撤退阻止 と、戦力の逐次投入でなくてもてる力をすべて出し切って日本を守ろうとした。 原発事故の最悪シナリオも専門家に依頼して作成して、それをパワーポイントでレクチャーさせて閣僚で危機感を共有していたという。 一方、安倍は菅直人が海水注入を止めたとデマを流しつづけた。名誉毀損裁判途中でデマメルマガをこっそりと削除。 今日も朝からクルクルパーのFラン統計ジジイが 他スレでお医者さんごっこだぁW 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 推定陽性者数: > summary(js$n.pos.sum) Min. 1st Qu. Median Mean 3rd Qu. Max. 3888 5094 5951 31450 8608 421156 > hdi(js$n.pos.sum)[1:2] lower upper 3888 172503 統計ジジイとコンプ薬屋はFラン仲間 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 >>109 あぁあ Fランジジイだぁ コテハンを命名してあげましょう コンプにならってコロナ薬屋とでも名乗るが良いですよ by 都内Sラン女子高生 今日も朝からクルクルパーのFラン統計ジジイが 他スレでお医者さんごっこだぁW 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 今日も朝からクルクルパーの Fラン統計ジジイ あらため コロナ薬屋が 他スレでお医者さんごっこだぁW 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 https://toyokeizai.net/sp/visual/tko/covid19/ のデータを使って P=14895/153581 # 2020/05/04 PCR検査の感度30-70%のモデルM1と感度70-90%のモデルM2のどちらが信憑性があるか、ベイズファクターで計算してみる。 M1は感度が最頻値60%標準偏差10%、M2は最頻値80%標準偏差10%のベータ分布に設定 特異度はいずれも最頻値95%標準偏差2.5%に設定し、有病率は一様分布を仮定 陽性数は、陽性率(P)=真陽性率+偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする。 事後確率分布は ://i.imgur.com/81bK4KE.png 陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと > d1/d2 # Savage-Dickey densiti ratio = BF12 [1] 1.007722 まあ、ちょっぴり、感度30-70%のモデルの方がいいかも、という結果。 陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。 味覚障害がある患者をPCR検査に回したら陰性であった。 PCR検査の感度は30-70%特異度は95-100%として、この患者が非感染者である確率の期待値と信頼区間を計算せよ。 (P.AA_COVID=1-72/202) # https://jamanetwork.com/journals/jama/fullarticle/2765183 sn=Mv2ab(0.5,0.2^2) sp=Mv2ab(0.95,0.05^2) data=list(sn=sn,sp=sp, P.AA_COVID=P.AA_COVID) modelString=' model{ # P(COVID | anosmia and/or ageusia) P.COVID_AA = P.AA_COVID*P.COVID/ (P.AA_COVID*P.COVID + P.AA_nonCOVID*(1-P.COVID)) # P(PCR- | aosmia and/or ageusia) P.nPCR = (1-P.COVID_AA)*spc + P.COVID_AA*(1-sen) # = TN + FN P.COVID_nPCR= (1-sen)*P.COVID/((1-sen)*P.COVID + spc*(1-P.COVID)) # P(COVID | PCR-) = P(PCR-|COVID)P(COVID)/P(PCR-) # = P(PCR-|COVID)P(COVID) / # (P(PCR-|COVID)P(COVID)+P(PCR-|nonCOVID)P(nonCOVID)) # prior P.AA_nonCOVID ~ dunif(0.056,0.127) # https://jamanetwork.com/journals/jama/fullarticle/2765183 P.COVID ~ dbeta(1,1) # prevalence sen ~ dbeta(sn[1],sn[2]) # sensitivity 0.3-0.7 spc ~ dbeta(sp[1],sp[2]) # specificity 0.9-1.0 } ' 今日も朝からクルクルパーの Fラン統計ジジイ あらため コロナ薬屋が 他スレでお医者さんごっこだぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 ボクちゃん ちゃんとコテハン付けるんでしよぉぉぉお 漢字はまだ読めないんでしたっけねえぇ by 都内Sラン女子高生 # s1,s2,n1,n2 -> BF10 BF01 <- function(s1,s2,n1,n2){ library(gmp) bf=chooseZ(n1,s1)*chooseZ(n2,s2)/chooseZ(n1+n2,s1+s2)*(n1+1)*(n2+1)/(n1+n2+1) BF01=as.numeric(bf) return(c('in favor of null-hyposisis : H0'=BF01)) } (bf01=BF01(s1=424,s2=5416,n1=777,n2=9072)) c('in favor of alternative hyposis : H1' = 1/as.numeric(bf01)) 今日も朝から人生が暇しかないの Fラン統計ジジイ あらため コロナ薬屋が 他スレでお医者さんごっこだぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 Fラン統計ジジイ あらため コロナ薬屋さん および コンプ薬屋さん って やっぱりバカなんですねぇぇえぇぇぇ! ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 Fラン統計ジジイ あらため コロナ薬屋さん 漢字は読めないんでしゅかぁぁあぁぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 ところで コンプ薬屋という方は 転売生活でもしてらっしゃるのでしょうか? http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/242 レスを遡ると品不足をアッピールして転売を煽るような書き込みをされてましたね ということは お仲間の Fラン統計ジジイあらため コロナ薬屋さんも ご職業は転売ヤーですか? by 都内Sラン女子高生 それで 転売ヤーの コロナ薬屋さん コンプ薬屋さんは 次はどんな品を転売予定ですか? それはそうと 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 library(rjags) s1=22;s2=10;n1=158;n2=78 c(s1/n1,s2/n2,s1/n1-s2/n2) shape1=shape2=0.5 data=list(s1=s1,s2=s2,n1=n1,n2=n2,shape1=shape1,shape2=shape2) modelString=' model{ # data s1 ~ dbin(p1,n1) s2 ~ dbin(p2,n2) # model p1 ~ dbeta(shape1,shape2) p2 ~ dbeta(shape1,shape2) delta = p1 - p2 # priors pri.delta= pri.p1 - pri.p2 pri.p1 ~ dbeta(shape1,shape2) pri.p2 ~ dbeta(shape1,shape2) } ' writeLines(modelString,'tmp.txt') n.iter=1e4 ; thin=1 jagsModel=jags.model('tmp.txt',data=data,n.chains=4,n.adapt=n.iter/5) update(jagsModel) codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin) plot(codaSamples) js=as.data.frame(as.matrix(codaSamples)) BEST::plotPost(js$delta,compVal = 0) library(polspline) fit=logspline(js$delta) pri.fit=logspline(js$pri.delta) xlim=c(-1,1) ; ylim=c(0,9) plot(fit,xlim=xlim, ylim=ylim,bty='l',lty=2) par(new=T) plot(pri.fit,xlim=xlim,ylim=ylim) dpost=dlogspline(0,fit) dprio=dlogspline(0,pri.fit) points(0,dpost) points(0,dprio,pch=19) dpost/dprio 転売ヤーの コロナ薬屋さん コンプ薬屋さん それが次の転売品ですか? 高音ぼったくりで通販サイトを荒らして 恥ずかしくないのですか? 転売品仕入れの時は夜中から並んで野糞をしてるって本当ですか? by 都内Sラン女子高生 http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/296 コロナ薬屋さん コンプ薬屋さんに共感されてますよ やっぱりお仲間だったんですね キンモー by 都内Sラン女子高生 今日も朝から人生が暇だけの コロナ薬屋 コンプ薬屋が 他スレでお医者さんごっこだぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 それで 私立薬学部卒の Fラン統計ジジイ あらため コロナ薬屋さん あなたと違って大学は旧帝受験予定ですが それはそうと昨夜は何を転売するため並んでたんですか? 野糞しながらぼったくり転売は恥ずかしくないのですか? by 都内Sラン女子高生 全てのパチンコ屋を営業停止にして コロナ患者の収容施設にしよう 3人の女子大生に18種類の体位を教えてtヶ月後に覚えていた体位の数は以下の通りである. ?は欠測値。 記憶保持率をexp(-α*t)+βのモデルに当て嵌めて 4人めの冬子を含めて欠測値を推定せよ。 t=1 2 4 7 12 21 35 59 99 200 治子 18 18 16 13 9 6 4 4 4 ? 奈津子 17 13 9 6 4 4 4 4 4 ? 亜希子 14 10 6 4 4 4 4 4 4 ? 冬子 ? ? ? ? ? ? ? ? ? ? それで 私立薬学部卒の Fラン統計ジジイ あらため コロナ薬屋さん お仕事は転売ヤーで間違いないですね 頭が悪いのに転売のお仕事の為に統計の勉強をされたんですね by 都内Sラン女子高生 統計上は PCR件数と 死亡数は 正の相関関係がある と なるほど PCR信者はバカなんですね by 都内Sラン女子高生 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 他スレでお医者さんごっこだぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 ネットを見ると、大学の9月入学を支持する高3生の泣き言、みたいなニュースがありますね 私の周りの高3生にはそんなことを支持する者は誰もいないのですが どんなFラン高の生徒の話なんでしょうね 受験は生モノ 勉強時間が足りないのや試験日に体調が悪くて不合格になるのは、本人の実力のうちでしょうに by 都内Sラン女子高生 えぇぇえぇぇぇっ Fラン統計ジジイ あらため コロナ薬屋さんって 反日パヨクなんですかぁぁあ? サイテーですね 下水のゾウリムシよりもはるかにサイッテーですね ダンゴムシのごとくゴロゴロ転がりながら 海より深く反省しなさい by 都内Sラン女子高生 >>130 わりと、面倒。 library(rjags) t <- c(1, 2, 4, 7, 12, 21, 35, 59, 99, 200) nt <- length(t) slist <- 1:4 ns <- length(slist) k <- matrix(c(18, 18, 16, 13, 9, 6, 4, 4, 4, NA, 17, 13, 9, 6, 4, 4, 4, 4, 4, NA, 14, 10, 6, 4, 4, 4, 4, 4, 4, NA, NA, NA, NA, NA,NA,NA,NA,NA,NA, NA), nrow=ns, ncol=nt, byrow=T) colnames(k)=paste0('t',as.character(t)) rownames(k)=paste0('subject',1:4) n <- 18 dataList=list(t=t,k=k,nt=nt,ns=ns,n=n) cat("model { for(i in 1:ns){ for(j in 1:nt){ k[i,j] ~ dbin(theta[i,j],n) predk[i,j] ~ dbin(theta[i,j],n) } } for(i in 1:ns){ for(j in 1:nt){ theta[i,j] <- min(1 - 1e-16, exp(-alpha[i]*t[j]) + beta[i]) } # min(1, ....) => Node inconsistent with parents -- ERROR } for(i in 1:ns){ alpha[i] ~ dbeta(shape1a,shape2a) # dnorm(mu.a,lambda.a)T(0,1) beta[i] ~ dbeta(shape1b,shape2b) # dnorm(mu.b,lambda.b)T(0,1) } shape1a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.a ~ dbeta(1,1) shape2a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.b ~ dbeta(1,1) shape1b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.a ~ dgamma(0.001,0.001)T(0.001,) shape2b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.b ~ dgamma(0.001,0.001)T(0.001,) } ", file='tmp.txt') pk4=NULL for(i in 1:10){ pk4=cbind(pk4,eval(str2lang(paste0("js$\'predk[4," ,i, "]\'")))) } colnames(pk4)=c(paste0('t0',1:9),'t10') boxplot(pk4) head(pk4) df4=as.data.frame(pk4) head(df4) predk4=tidyr::pivot_longer(df4,col=1:10) head(predk4) colnames(predk4)=c('time','items') library(ggplot2) g <- ggplot(predk4,aes(time,items)) g + geom_violin() g + geom_boxplot() reshape2 が tidyr::pivot_longer(df4,col=1:10) になっていた。 使っていないと、ggplot2も忘れてしまうなぁ。 g <- ggplot(predk4,aes(time,items)) g + geom_boxplot() https://i.imgur.com/krliw0d.png バイオリンプロットというよりナメクジプロットだな。 g + geom_violin() https://i.imgur.com/fYZHA53.png histの$densityを使って面積比で反映させてみた https://i.imgur.com/bIigw6R.png 標準のboxplotでも十分な気がする。 https://i.imgur.com/L210Sf8.png えぇぇえぇぇぇっ Fラン統計ジジイ あらため コロナ薬屋さん および コンプ薬屋さんって 裏口容疑者なんですかぁぁあ? だと思いましたぁぁあぁぁ 早く出頭するんですよ あなたたちのレスの全てから 頭の悪さを感じますからねえぇぇぇ by 都内Sラン女子高生 交絡因子を考慮しての記憶減衰の階層ベイズモデル デバッグできたので、事前分布をいろいろ変えてみて再現性を確認。 内視鏡休診で纏まった勉強する時間が取れて( ・∀・)イイ!! ##### psychophysical function with confounder cat(' model{ for(i in 1:nsubjs){ for(j in 1:nstim[i]){ z[i,j] ~ dbern(phi[i]) pi[i,j] ~ dunif(0,1) r[i,j] ~ dbin(theta[i,j],n[i,j]) theta[i,j] <- ifelse(z[i,j]==1,ilogit(alpha[i] + beta[i]*(x[i,j]-xmean[i])),pi[i,j]) } } for(i in 1:nsubjs){ phi[i] <- phi(probit) alpha[i] ~ dnorm(mua,lambdaa) beta[i] ~ dnorm(mub,lambdab) } probit ~ dnorm(muphi,lambdaphi) muphi ~ dnorm(0,0.001) mua ~ dnorm(0,0.001) mub ~ dnorm(0,0.001) lambdaphi ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0.3) lambdaa ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0,1000) lambdab ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0,1000) } ',file='tmp.txt') 開業医ではないあなたが 先に開業医板を荒らすから 私がここに降臨して差し上げてるんですよ それすら判らないくらい知能が低いんですかねぇ 女子高生にブチ切れる異常者のジジイって どんなホラー映画の地獄絵図ですかW by 都内Sラン女子高生 >>143 half cauchy にしても一様分布にしても事後分布はほとんど変化はないな。 モデルが優れているからだろう。 さて今日も学校は終わり。 お家に帰ってお勉強しましょう。 あ〜あ、Fラン統計ジジイは やっぱりくるくるぱーだったんですねぇ by 都内Sラン女子高生 # 相関係数の差の検定 n1=150 r1=0.29 n2=120 r2=0.5 pnorm(-abs(atanh(r1)-atanh(r2))/sqrt(1/(n1-3)+1/(n2-3))) 離散量のプロットは重なってしまうと相関がわかりにくくなる。 jitterつけると重なりは防げるな。 頻度に応じて面積比でプロット https://i.imgur.com/QVrj8hW.png n=10 ; r=500 x=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5)) y=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5)) plot(rep(0,24),pch=1:24,cex=2,axes=F,ann=F,ylim=c(0,0.1)) ;axis(1,at=1:24) xy=paste0(x,y) cex=numeric() for(i in 1:length(x)){ cex[i]=sum(xy[i]==xy) } layout(matrix(1:4,2,b=T)) par(mar=c(1,1,1,1)) plot(x,y,pch=16) plot(x,jitter(y),pch=16) plot(jitter(x),jitter(y),pch=16) plot(x,y,pch=16,cex=sqrt(cex)) とりあえず、完成。多分、もっと便利なパッケージがあるんだろうな。 graphics.off() oldpar=par() par(mar=c(1,1,1,1),oma=c(3,3,3,3),bty='l',ann=F,pch=16) plot(rep(0,24),pch=1:24,cex=2,axes=F,ann=F,ylim=c(0,0.1)) ;axis(1,at=1:24) n=10 ; r=250 x=sample(n,r,rep=T,prob=abs(rnorm(n))) y=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5)) plot(x,y) xy=paste0(x,y) cex=numeric() for(i in 1:length(x)){ cex[i]=sum(xy[i]==xy) } layout(matrix(1:4,2,b=T)) plot(jitter(x),jitter(y)) plot(x,y,cex=sqrt(cex)) plot(x,y,cex=3,col=rgb(0.1,0.1,0.1,0.25)) plot(x,y,cex=sqrt(cex),col=rgb(0.1,0.1,0.1,0.5)) layout(1) par(oldpar) 学校の先生に習いましたが、統計的にはPCR検査はむやみにやらないのが正解なんですね by 都内Sラン女子高生 【国民が誤解】安倍晋三「相談・受診の目安=PCR実施の基準ではない。国民への周知が足りなかった」 [スタス★] https://asahi.5ch.net/test/read.cgi/newsplus/1589284722/ 世界の頭脳は実態把握には大量検査が必要とわかっているなぁ。 東大の児玉名誉教授も同じ意見だな。 https://www.covid19-yamanaka.com/cont4/14.html 検診がメインだったから、内視鏡バイトは休診。 俺は6割の休業補償だが、銭ゲバ開業医はパート医を首にしているんだろうなぁ。 今日も朝から人生が暇しかない Fラン統計ジジイ あらため コロナ薬屋 が 発狂して連投だぁぁぁW ちゃんとコテハン付けないとダメでちよぉぉぉ 今思うのは反日野党が政権をとってたら 日本はとっくに壊滅していたであろうこと だって反日なんだもの by 都内Sラン女子高生 >>156 内視鏡バイトは休診なので、6割の休業補償はあるけど。 あんたのところは、感染拡大阻止に休診にして職員には休業補償してんの? まとまった時間がとれるので、https://bayesmodels.com/ の英文テキストで独学中。 行き詰まったら数学版で相談。 相関係数->FisherのZ変換->ベイズファクター r=js$r #The Fisher z-transformation: r.F <- atanh(r) fit.posterior <- logspline(r.F) plot(fit.posterior) posterior <- dlogspline(atanh(0), fit.posterior) r.F.prior <- atanh(runif(n.iter,-1,1)) fit.prior <- logspline(r.F.prior) prior <- dlogspline(atanh(0), fit.prior) BF10 <- prior/posterior BF10 まあ、ベイズファクターの計算結果はさほど変わらんな。 > BEST::plotPost(js$r,compVal = 0,showCurve = F,col='lightgreen') > lines(density(js$r),col='green') > abline(h=0.5,col=4) > fit.post=logspline(js$r) > points(0,dlogspline(0,fit.post),pch=19,col='white',cex=1.5) > points(0,dlogspline(0,fit.post),cex=1.5) > points(0,dunif(0,-1,1),pch=19,cex=1.5) > curve(dlogspline(x,fit.post),add=T) > fit.F=logspline(atanh(r)) > curve(dlogspline(x,fit.F), col=2, add=T) > d.prio=dunif(0,-1,1) > d.post=dlogspline(0,fit.post) > d.post/d.prio [1] 2.917773 > dlogspline(atanh(0),fit.F)/d.prio [1] 2.911407 あんまり相関係数の信頼区間など気にしていなかった。 必要ならmcmcかbootstrapすれば出せるし、 corr.testのソースを sink('cor.text.txt') stats:::cor.test.default sink() で覗いてみた。 if (n > 3) { if (!missing(conf.level) && (length(conf.level) != 1 || !is.finite(conf.level) || conf.level < 0 || conf.level > 1)) stop("'conf.level' must be a single number between 0 and 1") conf.int <- TRUE z <- atanh(r) sigma <- 1/sqrt(n - 3) cint <- switch(alternative, less = c(-Inf, z + sigma * qnorm(conf.level)), greater = c(z - sigma * qnorm(conf.level), Inf), two.sided = z + c(-1, 1) * sigma * qnorm((1 + conf.level)/2)) cint <- tanh(cint) FisherのZ変換して正規分布近似して信頼区間を出して、それをZ変換の逆関数で出しているな。 cat(' "名","本人","家族含む" "鳩山由紀夫首相",144269,144269 "菅直人副総理兼国家戦略担当相",905,2232 "原口一博総務相",914,1220 "千葉景子法相",3523,3523 "岡田克也外相",3273,8641 "藤井裕久財務相",14356,20214 "川端達夫文部科学相",4024,5583 "長妻昭厚生労働相",0,891 "赤松広隆農相",4864,5934 "直嶋正行経済産業相",3333,3333 "前原誠司国土交通相",741,1441 "小沢鋭仁環境相",2089,4014 "北沢俊美防衛相",309,609 "平野博文官房長官",1195,1875 "中井洽国家公安委員長",1296,1296 "亀井静香金融・郵政担当相",9427,18745 "福島瑞穂消費者・少子化担当相",12734,25000 "仙谷由人行政刷新担当相",1968,3987 ', file='shisan.csv') X=read.csv("shisan.csv") X library(boot) md = boot(X, function(df,idx) median(df[idx,2]), R=1e4) plot(md) mn = boot(X, function(df,idx) mean(df[idx,2]), R=1e4) plot(md) boxplot(list(median=md$t,mean=mn$t),horizontal = TRUE, las=1) >>162 パッケージ bootをつかうと library(boot) df=data.frame(x,k) bt=boot(df,function(df,idx) cor(df[idx,1],df[idx,2]),R=1e5) plot(bt) HDInterval::hdi(bt$t)[1:2] BEST::plotPost(as.vector(bt$t),xlab=bquote(cor),compVal = 0) > HDInterval::hdi(bt$t)[1:2] [1] -0.06736828 0.31237556 https://i.imgur.com/vCeo1Yn.png >>156 アベノミクスでの経済成長率は震災があった民主の頃より低いという現実を認められないんだなぁ。 日本を破壊したいなら、安倍一択。 安倍の長期政権で中国も韓国も笑いが止まらないはず、 日本が衰退する政策ばかりを選択しているから。 野党が反感をもたれて安倍信者が増えることこそ日本破壊が進んで パヨクの思う壺 結局、これが核心部分 Rt =I t (number of new infections generated at time t) / Σ[s=1,t] I t-s * Ws ( = total infectiousness of infected individuals at time t) Ws : an infectivity profile given by a probability distribution ws, dependent on time since infection of the case, s, but independent of calendar time, t. E[I t] = Rt Σ[s=1,t] I t-s * Ws Σ[s=1,t] I t-s * Wsの部分が畳み込み積分で Ws ∝ serial interval # https://www.nejm.org/doi/full/10.1056/NEJMoa2001316 #--- incubation period --- # from Li et al NEJM 2020 # lognormal mean = 5.2 ln.par1 = 1.434065 ln.par2 = 0.6612 x=rlnorm(1e6,ln.par1,ln.par2) BEST::plotPost(x) curve(dlnorm(x,ln.par1,ln.par2),add=T) quantile(x,c(0.025,0.5,0.975)) " c.onset = c.infected + c.incubation d.onset = d.infected + d.incubation c.onset - d.onset = onset.delay c.infected - d.infected + c.incubation - d.incubation = delay c.infected - d.infected = onset.delay + d.incubation - c.incubation " onset.delay = 2 d.incubation = rlnorm(1e6,ln.par1,ln.par2) c.incubation = rlnorm(1e6,ln.par1,ln.par2) infection.delay = onset.delay + d.incubation - c.incubation BEST::plotPost(infection.delay,compVal = 0) read.cgi ver 07.5.5 2024/06/08 Walang Kapalit ★ | Donguri System Team 5ちゃんねる