判別分析の判別線がおかしい

という相談を受けた。
解析したいデータは2変数なので、例として
判別分析モデルの応用
がついてきたが、3変数の場合はirisを用いて
データ解析・マイニングとR言語
もついてきて、こんな感じでやりたい、と言われた。

0/1のtype は以下の線形モデルで表現されて
Type = \alpha + \beta_xX+\beta_yY
判別線は
Y=-(\alpha + \beta_xX)/\beta_y
となる。

しかし、リンクの通りに\alphaを計算すると、判別線より上か下かを決めるスコア(正負で上下が決まる)と、判別線との関係がなにかおかしい、というのが問題のようである。
図で言うと、青の判別線の下にある8個の☓がおかしい点である。
f:id:MikuHatsune:20210212164704p:plain

結論から言うと

# 以上の手順は1回で済ますこともできる
apply(disc.model$means %*% disc.model$scaling, 2, mean)

の部分が平均を取っており、リンクの例題なら2カテゴリーが1:1 で作成されたので偶然結果が一致するが、今回のデータは0/1 が偏っているので、その偏りに応じた重み付けが必要になる。
これはiris を使った2番目のリンクの例題でも、1:1:1 になって平均を取って偶然一致するので、わかりにくい。

alpha <- c(disc.model$prior %*% disc.model$means %*% disc.model$scaling)

とすればデータの事前確率に応じた重み付けが取れる(行列演算なので返り値がベクトルになるようにcをかませている)。

緑の判別線ならうまくいったようである。
f:id:MikuHatsune:20210212165254p:plain

txt <- "
type X Y
0 0 107
0 0 109
0 0 97
0 0 75
0 0 101
0 0 92
0 0 88
0 0 99
0 0 89
0 22 96
0 37 110
0 0 90
0 0 101
0 92 68
0 51 99
0 26 105
0 0 92
0 0 81
0 0 95
0 0 96
0 18 88
0 0 89
0 0 94
0 0 92
0 0 90
0 0 109
0 0 88
0 0 84
0 0 108
0 0 92
1 0 75
1 47 61
1 29 73
1 96 76
1 46 85
1 103 72
1 68 60
1 0 78
1 0 52
1 0 84
1 0 64
1 64 85
1 72 89
1 50 65
1 50 75
1 116 87
1 79 80
1 130 77
1 48 82
1 40 86
1 140 71
1 34 77
1 109 95
1 208 67
1 85 81
1 0 62
1 91 64
1 106 95
1 27 73
1 36 61
1 28 77
1 53 94
1 113 94
1 59 77
1 38 63
1 0 75
1 78 60
1 321 68
1 0 55
1 81 54
1 77 87
1 49 57
1 131 69
1 100 54
1 149 51
1 59 65
1 83 43
1 39 56
1 176 47
1 122 60
1 172 78
1 123 46
1 164 54
1 42 64
1 86 60
1 31 75
1 27 73
1 64 57
1 124 72
1 59 62
1 0 84
1 18 66
1 103 92
1 0 76
1 22 44
1 125 54
1 197 71
"

dat <- read.table(text=txt, header=TRUE)

library(MASS)
disc.model <- lda(type ~ Y + X, data=dat)
p <- predict(disc.model)
pp <- p$posterior
p2 <- function(x) ifelse(x < 1, 1/x, x)

TF <- p$class == dat$type

group.mean <- disc.model$mean
# 線形係数のデータ行列をcoefに代入
coef <- disc.model$scaling
# 行列のかけ算は%*%で行える
# その結果をaに代入
a <- group.mean %*% coef
# 行列aの列平均を計算する
# 引数に2ではなく1を指定すれば行に対して計算される
# 仮にmeanをsumにすると、平均ではなく合計が計算される
apply(a, 2, mean)
# 以上の手順は1回で済ますこともできる
# web のやつ
alpha1 <- apply(disc.model$means %*% disc.model$scaling, 2, mean)
# データの事前確率を反映させた正しいもの
alpha2 <- c(disc.model$prior %*% disc.model$means %*% disc.model$scaling)

# 判別線
f2 <- function(x, alpha) -(-alpha + disc.model$scaling["X", "LD1"]*x)/disc.model$scaling["Y", "LD1"]
x <- seq(0, 500, length=1000)

plot(Y ~ X, data=dat, col=dat$type+1, pch=ifelse(p$x < 0, 4, 16))
legend("topright", legend=sprintf("type %d", 0:1), pch=16, col=1:2)
lines(x, f2(x, alpha1), lty=3, col=4, lwd=3)


plot(Y ~ X, data=dat, col=dat$type+1, pch=ifelse(p$x < 0, 4, 16))
lines(x, f2(x, alpha1), lty=3, col=4, lwd=3)
lines(x, f2(x, alpha2), lty=3, col=3, lwd=3)
legend("topright", legend=sprintf("type %d", 0:1), pch=16, col=1:2)
legend("bottomright", legend=c("正", "負"), pch=c(4, 16), col=1:2)


plot(Y ~ X, data=dat, col=dat$type+1, pch=c(4, 16)[TF+1], cex=log10(p2(pp[,1]/pp[,2]))+1)
legend("topright", legend=sprintf("type %d", 0:1), pch=16, col=1:2)
lines(x, f2(x, alpha2), lty=3, col=4, lwd=3)

新型肺炎の年齢階級別死亡率

新型コロナウイルス 国内感染の状況
こちらの2021年1月13日時点の各年代別感染者数と死亡者数のデータがあったので単純に推定してみる。

その1
drcパッケージにあるLL.4関数は、4パラメータモデルのロジスティック曲線に当てはめる。
これで簡単にやる

その2
rstan を使う。年齢を経るごとに死亡率が上がるとして、orderedで昇順の値のサンプリングをして、inv_logitにより0~1 の値になるようにする。すなわちこれが死亡率になり、やっていることはロジスティック曲線への変換である。

死亡率は以下のようになる(%)

      age             
 [1,]   5  0.007447185
 [2,]  15  0.007450361
 [3,]  25  0.007735825
 [4,]  35  0.013074665
 [5,]  45  0.059155769
 [6,]  55  0.310493872
 [7,]  65  1.314663597
 [8,]  75  4.443249524
 [9,]  85 12.025304864

f:id:MikuHatsune:20210118232728p:plain

rstan でやるよりかは、推定曲線へのフィッティングのほうが、predict関数でinterval="confidence"を指定すると区間推定値が出るので有用だと思う。
例えば100歳での死亡率は95%信頼区間

     Lower Prediction      Upper 
  30.62692   33.85670   37.08648 

となる。

# 1/13
death <- cbind(c(0, 0, 2, 11, 35, 113, 328, 944, 2401),
               c(7160, 18377, 69311, 46455, 42939, 39148, 24823, 21243, 19967))
rownames(death) <- c("10代以下", sprintf("%d代", seq(10, 70, by=10)), "80代以上")

age <- seq(5, 85, by=10)
p <- death[,1]/death[,2]



library(vioplot)
library(drc)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

code <- "
data{
  int N;
  int death[N];
  int positive[N];
}
parameters{
  ordered[N] theta;
}
model{
  for(i in 1:N){
    death[i] ~ binomial(positive[i], inv_logit(theta[i]));
  }
}
"

m0 <- stan_model(model_code=code)
standata <- list(N=nrow(death), death=death[,1], positive=death[,2])
fit <- sampling(m0, data=standata, iter=3000, warmup=1500, seed=1234)

ex <- extract(fit)
r <- 1/(1+exp(-ex$theta))

d <- drm(p ~ age, fct=LL.4())
x <- seq(0, 100, by=1)
y <- predict(d, as.data.frame(x), interval="confidence")

xl <- range(x)
yl <- range(0, r, y)

plot(d, log="", xlim=xl, ylim=yl, xlab="年齢", ylab="陽性者死亡率", bty="n", axes=FALSE)
lines(x, y[,2], col=2)
lines(x, y[,3], col=2)
vioplot(r, add=TRUE, at=age, colMed="yellow", wex=5)
axis(1, lwd=0, lwd.ticks=1, at=age, labels=age)
axis(2, lwd=0, lwd.ticks=1)

高水準関数の作図後に低水準関数で作図したときの線の扱いをなんとかしたい

意味的にはこんな感じである。
plotで作図すると高水準関数といって、とりあえず作図されるが、その後に低水準関数と呼ばれる、例えばpointslineといった関数で後付でお絵描きができる。
自分がよくやるのはplotしたあとに、ある領域はこうである、みたいな色付けをpolygon関数でやるのだが、何も考えないでplotのあとにpolygonを行うとこんな感じになる。

d <- 0.5
plot(1, xlab="", ylab="")
box(lwd=10)
pa <- par()$usr
polygon(c(pa[1]+d, pa[2]-d, pa[2]-d, pa[1]+d), rep(pa[3:4], each=2), border=NA, col=grey(0.8))

f:id:MikuHatsune:20210107143612p:plain
box関数で枠の太さを水増ししているが、要はpolygonで作図した線がplotで定義されている枠線のちょうど真ん中を通るように設計されているので、polygonがレイヤーでいうところの最上位に来ていてこんな感じのなんか気持ち悪い感じになる。

これを解決してpolygonが下に来てplotの枠線が上に来るようにするには、それはそうするしかない。が、普通にboxで後付けするとplotで描かれる枠線とかぶって更に気持ち悪くなる。
というわけでplotするときに、不要なものやあとで付け足したいものは引数で描かないように設定する。

plot(1, xlab="", ylab="", frame=FALSE, xaxt="n", yaxt="n")
axis(1, lwd=0, lwd.ticks=1)
axis(2, lwd=0, lwd.ticks=1)
polygon(c(pa[1]+d, pa[2]-d, pa[2]-d, pa[1]+d), rep(pa[3:4], each=2), border=NA, col=grey(0.8))
box(lwd=10)

f:id:MikuHatsune:20210107143950p:plain

メタアナリシスっぽいので公平な入試を受けたい 令和2年版

この記事はR Advent Calendar 2020 - Qiitaの24日の配当記事です(書くネタなかったけど無理やりこじつけた
メタアナリシスの解析にmeta パッケージを使っています。


令和2年版の調査結果というものを観測した。
医学部医学科の入学者選抜における公正確保等に係る調査について:文部科学省

いままでと同じようにやってみる。
mikuhatsune.hatenadiary.com
mikuhatsune.hatenadiary.com

結論から言うと男性が1.1倍程度合格しやすいようである。
f:id:MikuHatsune:20201226002714p:plain

library(meta)
dat <- read.csv(text="
大学,男性受験者数,男性合格者数,女性受験者数,女性合格者数
北海道大学, 266,  89, 77, 25
旭川医科大学, 336, 224, 66, 37
弘前大学, 225, 155, 66, 46
東北大学, 294, 109, 89, 32
秋田大学, 253, 185, 82, 48
山形大学, 212, 179, 63, 44
筑波大学, 238, 208, 91, 49
群馬大学, 205, 130, 82, 35
千葉大学, 260, 116, 90, 39
東京大学, 271,  67, 81, 20
東京医科歯科大学, 244, 130, 73, 43
新潟大学, 297, 121, 93, 30
富山大学, 218, 165, 54, 52
金沢大学, 204,  70, 89, 24
福井大学, 205, 145, 75, 40
山梨大学, 274,  88,111, 32
信州大学, 219, 164, 71, 55
岐阜大学, 369, 149, 76, 29
浜松医科大学, 213, 141, 78, 42
名古屋大学, 206, 101, 79, 32
三重大学, 213, 124, 83, 42
滋賀医科大学, 169, 133, 45, 51
京都大学, 215,  64, 91, 17
大阪大学, 235,  76, 84, 17
神戸大学, 164,  96, 71, 41
鳥取大学, 398, 250, 66, 39
島根大学, 272, 238, 59, 45
岡山大学, 242, 169, 65, 46
広島大学, 335, 154, 81, 39
山口大学, 295, 155, 70, 38
徳島大学, 125,  75, 73, 45
香川大学, 264, 219, 69, 45
愛媛大学, 318, 277, 57, 58
高知大学, 271, 238, 69, 44
九州大学, 200,  55, 91, 20
佐賀大学, 229, 149, 60, 43
長崎大学, 223, 125, 85, 36
熊本大学, 358, 179, 79, 33
大分大学, 255, 196, 61, 41
宮崎大学, 179, 147, 69, 45
鹿児島大学, 233, 154, 72, 40
琉球大学, 229, 156, 64, 50
札幌医科大学, 249, 124, 81, 31
福島県立医科大学, 207, 138, 89, 48
横浜市立大学, 130, 110, 58, 42
名古屋市立大学, 135,  98, 61, 37
京都府立医科大学, 158,  80, 79, 29
大阪市立大学, 145,  90, 70, 25
奈良県立医科大学, 435, 264, 84, 40
和歌山県立医科大学, 131,  47, 78, 23
岩手医科大学,1620, 836,190, 72
東北医科薬科大学,1060, 522,228, 97
自治医科大学,1604,1030, 75, 48
獨協医科大学,2081,1310,179, 94
埼玉医科大学,2753,1972,144,105
国際医療福祉大学,2018,1769,230,161
杏林大学,1916,1267,110, 81
慶應義塾大学, 841, 332,123, 44
順天堂大学,2164,1607,200,125
昭和大学,2438,1661,166,119
帝京大学,5263,3607,182,125
東京医科大学,1454, 923,183,107
東京慈恵会医科大学,1251, 549,288,114
東京女子医科大学,NA,1478,NA,343
東邦大学,1484,1121, 69, 61
日本大学,1784,1160,128, 66
日本医科大学,2250,1470,324,220
北里大学,1241, 755,182, 97
聖マリアンナ医科大学,1169, 965, 71, 96
東海大学,2608,1838,146, 87
金沢医科大学,2789,1725,128,109
愛知医科大学,2084,1370,226,130
藤田医科大学,2041,1233,210,120
大阪医科大学,1681,1043,126, 86
関西医科大学,2239,1591,354,223
近畿大学,2060,1353,159, 85
兵庫医科大学,1155, 912,137,123
川崎医科大学,1034, 623, 71, 55
久留米大学,1949,1100,127, 71
産業医科大学, 928, 621, 76, 46
福岡大学,2233,1363,110, 65"
)

b <- metabin(女性受験者数, 男性受験者数, 女性合格者数, 男性合格者数, data=dat, sm="OR", studlab=大学, comb.fixed=TRUE, comb.random=TRUE, label.e="男性", label.c="女性")

alpha <- 0.05
cols <- ifelse(b$pval < alpha, ifelse(b$TE < 0, "red", "blue"), "black")

forest(b, label.right="男性有利", label.left="女性有利",
       col.label.right="blue", col.label.left="red",
       col.study=cols, lwd=2)

虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 12話

侑と歩夢回だったが孤立したノードがたくさん。
f:id:MikuHatsune:20201220203756p:plain

虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 1〜3話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 4話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 5話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 6話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 7話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 8話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 9話と10話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 11話 - 驚異のアニヲタ社会復帰の予備

虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 11話

歩夢のレズメンヘラっぷりが炸裂したが果林としずくが相変わらずハブられていた。
f:id:MikuHatsune:20201213144205p:plain
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 1〜3話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 4話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 5話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 6話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 7話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 8話 - 驚異のアニヲタ社会復帰の予備
虹ヶ咲学園スクールアイドル同好会の名前呼び合いグラフをかく 9話と10話 - 驚異のアニヲタ社会復帰の予備