現在のアクセス件数は24万3千件です。
必要としていただいている皆様のために、只今リニューアルの準備中です。
しばらく（たぶん半年くらい・・・）投稿や更新が滞りますので予めご了承ください。
中心極限定理
投稿日:2016/10/31, 最終更新日:2022/3/23
・母集団の分布がどのような分布であっても、無作為抽出した標本における和の分布は、標本の大きさｎが大きいときに正規分布に収束する。
= 母集団の分布がどのような分布であっても、無作為抽出した標本における標本平均の分布は、標本の大きさｎが大きいときに正規分布に収束する。
表の確率1/3のコインでの実験
このコインを数千回投げる実験です。投げる枚数を徐々にふやして表が出る枚数の平均値を求めます。でも実際に投げるのは大変なので、2項分布を利用して数千回投げた場合の表の出る確率を求めます。
実験1) 1枚のコインを数千回投げた場合に表が0枚, 1枚となる確率
x<-c() s<-0 q<-1 for(i in s:q)x <-c (x, i) #空のベクトルx（）にiをくっつける ベクトルX=c(0, 1)の完成 # x # 0 1
確率 に従う二項分布より確率を求める
n<-1 #コインを１枚投げた場合 p<-1/3 # 表の出る確率 x1 <- choose(1, x)*(p^x)*(1 - p)^(n - x) #表が0枚, 1枚となる確率 #累積確率 cumsum(x1) #グラフを描く barplot(x1, space=c(0.5), names.arg=c("0枚", "1枚"))
実験2) 2枚のコインを数千回投げた場合に表が0枚, 1枚, 2枚 となる確率
x <- c() s <- 0 q <- 2 for(i in s:q)x <- c(x, i) x2 <- choose(2, x)*(p^x)*(1 - p)^(n - x) cumsum(x2) barplot(x2,space=c(0.5),names.arg=c("0枚", "1枚","2枚")) #二項分布を示している
実験3) 5枚のコインを数千回投げた場合に表が0･･･5枚 となる確率
x <- c() s <- 0 q <- 5 for(i in s:q)x <- c(x, i) x3 <- choose(5, x)*(p^x)*(1 - p)^(n - x) cumsum(x3) barplot(x3,space=c(0.5),names.arg=c("0枚", "1枚","2枚","3枚","4枚","5枚"))
実験4) 10枚のコインを数千回投げた場合に表が0･･･10枚 となる確率
x <- c() s <- 0 q <- 10 for(i in s:q)x <- c(x, i) x4 <- choose(10, x)*(p^x)*(1 - p)^(n - x) cumsum(x4) barplot(x4, space=c(0.5), names.arg=c(0:10)) #「表が出る枚数の平均値」が正規分布に10/3 の付近に収束していることが理解できる
実験5) 100枚のコインを数千回投げた場合に表が0･･･100枚 となる確率
x <- c() s <- 0 q <- 100 for(i in s:q)x <- c(x, i) x5 <- choose(100, x)*(p^x)*(1 - p)^(n - x) cumsum(x5) barplot(x5, space=c(0.5), names.arg=c(0: 100)) #「表が出る枚数の平均値」が100/3に収束して、正規分布を示している。
実験で描いたグラフ
par(mfrow = c(2,3)) barplot(x1, space=c(0.5), names.arg = c("0枚", "1枚"), main="実験1") barplot(x2, space=c(0.5), names.arg = c("0枚", "1枚","2枚") , main="実験2") barplot(x3, space=c(0.5), names.arg = c("0枚", "1枚","2枚","3枚","4枚","5枚"), main="実験3") barplot(x4, space=c(0.5), names.arg = c(0:10), main="実験4") barplot(x5, space=c(0.5), names.arg = c(0:100), main="実験5") par(mfrow = c(1, 1))
Nagelkerke's R squared
Nagelkerkes's R2 の求め方
logL: 最大対数尤度
deviance ＝-2*logL
Dnull: Null modelのdeviance
D: 予測したモデル（変数あり）のdeviance
n: サンプルサイズ
R2 : Nagelkerkes's R2
irisを使用してRのパッケージ(fmsb)で算出
library(fmsb) dat <- iris[iris$Species=="setosa" | iris$Species=="virginica" ,] dat$Species <- ifelse(dat$Species=="setosa",1,0) fit_null <- glm(Species ~ 1, family = binomial (link = "logit"), data = dat) fit <- glm(Species ~ Sepal.Width, family = binomial (link = "logit"), data = dat) NagelkerkeR2(fit) # N=100 # R^2 = 0.4017335
求め方
n = dim(dat)[1] #最大対数尤度 L0 <- logLik(fit_null)[1] L1 <- logLik(fit)[1] x <- 1- (exp(L0)/exp(L1))^(2/n) y <- 1- (exp(L0))^(2/n) x/y # R^2 = 0.4017335 #Rを使ってdevianceから求めると・・・ Dnull = summary(fit)$null.deviance D = summary(fit)$deviance x <- 1 - exp((D - Dnull)/n) y <- 1 - exp(-Dnull/n) x/y # R^2 = 0.4017335
参考
* E.W.SteyerbergClinical; Prediction Models, 2nd ed , A Practical Approach to Development, Validation, &Updating, 2019, p66
* 久保 拓弥; データ解析のための統計モデリング入門―一般化線形モデル・階層ベイズモデル・ＭＣＭＣ, 2012, pp94-111
* FAQ: What are pseudo R-squareds?
最小二乗法と最尤推定法
投稿日:2022.3.9
irisを使用した単回帰分析
まずは最小二乗法
dat = iris fit1 <- lm(Sepal.Length ~ Sepal.Width, data = dat) summary(fit1) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5262 0.4789 13.63 <2e-16 *** Sepal.Width -0.2234 0.1551 -1.44 0.152 --- Residual standard error: 0.8251 on 148 degrees of freedom Multiple R-squared: 0.01382, Adjusted R-squared: 0.007159 F-statistic: 2.074 on 1 and 148 DF, p-value: 0.1519
残差の標準誤差や決定係数などが算出されます。
またF値やp値は次に示すような分散分析の結果になります。
anova(fit1) Df Sum Sq Mean Sq F value Pr(>F) Sepal.Width 1 1.412 1.41224 2.0744 0.1519 Residuals 148 100.756 0.68078
参考
次に最尤推定法による単回帰分析（lmがglmに変わるだけ・・・です）。もちろん答えは同じです。
fit2 <- glm(Sepal.Length ~ Sepal.Width, data = dat) summary(fit2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5262 0.4789 13.63 <2e-16 *** Sepal.Width -0.2234 0.1551 -1.44 0.152 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 0.6807844) Null deviance: 102.17 on 149 degrees of freedom Residual deviance: 100.76 on 148 degrees of freedom AIC: 371.99
ただし対数尤度を利用した解析結果ですのでdeviance、AICが算出されます。
確率分布は正規分布で、リンク関数は恒等リンクがデフォルトなので、本当は以下のように書かれています（算出結果は全く同じ）。
fit3 <- glm(Sepal.Length ~ Sepal.Width, family = gaussian (link = identity), data = dat) summary(fit3)
family = gaussianのリンク関数は、 identity, log, inverseが使用できるようです。
The gaussian family accepts the links (as names) identity, log and inverse
ちなみにリンク関数を変えたら、（当たり前ですが）異なる結果となります。
fit4 <- glm(Sepal.Length ~ Sepal.Width, family = gaussian(link=log), data = dat) summary(fit4) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.87900 0.08200 22.914 <2e-16 *** Sepal.Width -0.03723 0.02668 -1.396 0.165 (Dispersion parameter for gaussian family taken to be 0.6810537) Null deviance: 102.17 on 149 degrees of freedom Residual deviance: 100.80 on 148 degrees of freedom AIC: 372.05
連続変数をダミー変数に変換する方法
投稿日：2019-03-13、最終更新日：2022.2.15
makedummiesの作者の荒様から,便利な使い方を教えていただきましたので再々更新しておきます。
irisのSepal.Lengthを使用して、4～5, 5～6, 6～7, 7～8 の4区分をダミー変数にします。
一瞬で終わります
library(dplyr) library(makedummies) data <- data.frame(iris$Sepal.Length) summary(data) iris.Sepal.Length Min. :4.300 1st Qu.:5.100 Median :5.800 Mean :5.843 3rd Qu.:6.400 Max. :7.900 # 4～8の間なので4以上5未満, 5以上6未満, 6以上7未満, 7以上8未満 の4区分でダミー変数を作成します dat <- iris %>% mutate( dummy = cut(Sepal.Length, breaks = c(4,5,6,7,8), right = FALSE), dummy = factor(dummy, labels =c("4-5", "5-6", "6-7", "7-8"))) #分かりにくいので1列,５列とdummy変数のみを残します dat <- dat[,c(-2,-3,-4),] Sepal.Length Species dummy 1 5.1 setosa 5 2 4.9 setosa 4 3 4.7 setosa 4 #このままダミー変数にしたら名義変数を全部ダミーにしてしまいます。なので、ダミー変数にしたくない名義変数はas.is = ""としておきます。 #基準を作りたくないとき（全部ダミー変数にしたいとき）には basal_level = T dat <- makedummies(dat, basal_level = F, as.is = "Species") Sepal.Length Species dummy_5-6 dummy_6-7 dummy_7-8 1 5.1 setosa 1 0 0 2 4.9 setosa 0 0 0 3 4.7 setosa 0 0 0 4 4.6 setosa 0 0 0 5 5.0 setosa 1 0 0 6 5.4 setosa 1 0 0 ＊4以上5未満が基準になっているので、変数としては出てきません
ここからは個人的な学習用です
備忘録として、パイプを分解して解説入れておきます
# dplyr::mutate データセットに新たに変数を追加する関数です #変数dummyを追加します dat2 <- mutate(iris, dummy = cut(Sepal.Length, breaks = c(4,5,6,7,8), right = FALSE)) head(dat2) Sepal.Length Sepal.Width Petal.Length Petal.Width Species dummy 1 5.1 3.5 1.4 0.2 setosa [5,6) 2 4.9 3.0 1.4 0.2 setosa [4,5) 3 4.7 3.2 1.3 0.2 setosa [4,5) 4 4.6 3.1 1.5 0.2 setosa [4,5) 5 5.0 3.6 1.4 0.2 setosa [5,6) 6 5.4 3.9 1.7 0.4 setosa [5,6) dat3 <- mutate(dat2, dummy = factor(dummy, labels =c("A", "B", "C", "D"))) #カテゴリーなので数字じゃなくても良い head(dat3) # 変数dummyをカテゴリー化 Sepal.Length Sepal.Width Petal.Length Petal.Width Species dummy 1 5.1 3.5 1.4 0.2 setosa B 2 4.9 3.0 1.4 0.2 setosa A 3 4.7 3.2 1.3 0.2 setosa A 4 4.6 3.1 1.5 0.2 setosa A 5 5.0 3.6 1.4 0.2 setosa B 6 5.4 3.9 1.7 0.4 setosa B # ダミー変数 dat4 <- makedummies(dat3, basal_level = TRUE, as.is = "Species") head(dat4) Sepal.Length Sepal.Width Petal.Length Petal.Width Species dummy_A dummy_B dummy_C dummy_D 1 5.1 3.5 1.4 0.2 setosa 0 1 0 0 2 4.9 3.0 1.4 0.2 setosa 1 0 0 0 3 4.7 3.2 1.3 0.2 setosa 1 0 0 0 4 4.6 3.1 1.5 0.2 setosa 1 0 0 0 5 5.0 3.6 1.4 0.2 setosa 0 1 0 0 6 5.4 3.9 1.7 0.4 setosa 0 1 0 0 # 基準の列を除いたダミー変数 dat4_2 <- makedummies(dat3, as.is = "Species") head(dat4_2) Sepal.Length Sepal.Width Petal.Length Petal.Width Species dummy_B dummy_C dummy_D 1 5.1 3.5 1.4 0.2 setosa 1 0 0 2 4.9 3.0 1.4 0.2 setosa 0 0 0 3 4.7 3.2 1.3 0.2 setosa 0 0 0 4 4.6 3.1 1.5 0.2 setosa 0 0 0 5 5.0 3.6 1.4 0.2 setosa 1 0 0 6 5.4 3.9 1.7 0.4 setosa 1 0 0
やっぱり便利なパッケージです
ここからは、以前の記事です
大変な作業でした。
適当に150人分の年齢データを作成します
乱数で作成しているので、ファイルとして保存します
何度読み込んでも同じデータで練習できます
as.integer(runif(50 ,min =30 ,max=51)) as.integer(runif(50 ,min =40 ,max=80)) as.integer(runif(50 ,min =60 ,max=90)) age <- c(as.integer(runif(50 ,min =30 ,max=51)), as.integer(runif(50 ,min =40 ,max=80)), as.integer(runif(50 ,min =60 ,max=90))) da <- data.frame(age)
write.csv(da, file="age.csv")
保存したファイルを読み込み込んで作業を始めます
data_age <- read.csv("age.csv")
以下の3群のダミー変数を作成してみます
30～49歳
50～69歳
70～80歳
#30歳～50歳に該当する行のみTRUEとなる a30 <- data_age$age <=49 #30歳～50歳に1、他は0とするダミー変数を作成する age_30c <- rep(1,length(data_age[a30,"age"])) age_50c <- rep(0,length(data_age[a30,"age"])) age_70c <- rep(0,length(data_age[a30,"age"])) #3列追加する 列名をage30d、age40d、age50dとする a30_data <- transform(data_age[a30,],age30d=age_30c, age50d=age_50c,age70d=age_70c ) a50 <- data_age$age >=50 & data_age$age <=69 age_30c <- rep(0,length(data_age[a50,"age"])) age_50c <- rep(1,length(data_age[a50,"age"])) age_70c <- rep(0,length(data_age[a50,"age"])) a50_data <- transform(data_age[a50,],age30d=age_30c, age50d=age_50c,age70d=age_70c ) a70 <- data_age$age >=70 age_30c <- rep(0,length(data_age[a70,"age"])) age_50c <- rep(0,length(data_age[a70,"age"])) age_70c <- rep(1,length(data_age[a70,"age"])) a70_data <- transform(data_age[a70,],age30d=age_30c, age50d=age_50c,age70d=age_70c ) #ダミー変数が追加されたデータセットをdata_age_dとして作成 data_age_d <- rbind(a30_data, a50_data, a70_data)
完成したデータセットにイメージ
data_age_d[0:10,] X age age30d age50d age70d 1 1 32 1 0 0 2 2 39 1 0 0 3 3 43 1 0 0 4 4 34 1 0 0 5 5 32 1 0 0 6 6 46 1 0 0 7 7 43 1 0 0 8 8 35 1 0 0 9 9 34 1 0 0 10 10 35 1 0 0
さらに簡単に入力
30歳～49歳を対象としたサンプルを想定
30代をreferenceとして解釈したい場合は…
dataset=dat , 変数名=年齢
y <- dat x <- dat$年齢 aa <- x<=39 #ここが基準になるので、全部0 a1 <- rep(0,length(y[aa, "年齢"])) AA <- transform(y[aa,], age30d=a1, age40d=a1, age50d = a1 ) bb <- x >=40 & x <=49 b1 <- rep(1,length(y[bb, "年齢"])) b2 <- rep(0,length(y[bb, "年齢"])) BB <- transform(y[bb,], age30d=b2, age40d=b1, age50d = b2 ) cc <- x >=50 c1 <- rep(1,length(y[cc, "年齢"])) c2 <- rep(0,length(y[cc, "年齢"])) CC <- transform(y[cc,], age30d=c2, age40d=c2, age50d = c1 ) data_age_d <- rbind(AA, BB, CC)
ヒストグラム
library(ggplot2)
data <- c(rnorm(n=100, mean=5, sd=1),rnorm(n=100, mean=3, sd=3)) cat <- c(rep("x",100), rep("y",100)) dat1 <- data.frame(cat, data) colnames(dat1) <- c("name","value") #乱数なのでデータは変わります
単純なヒストグラム
HG1= ggplot(dat1, aes(x = value)) + geom_histogram(binwidth = 1) plot(HG1)
color=枠、fill=塗りつぶし、alpha=透かし (0-1)
HG2= ggplot(dat1, aes(x = value)) + geom_histogram(binwidth = 1, color="black", fill="grey", alpha=0.6) plot( HG2 )
各グループで色分け、凡例の位置（position）
HG3= ggplot(dat1, aes(x = value, color = name, fill =name)) + geom_histogram(position = "identity", binwidth = 1, alpha=0.6) plot(HG3)
枠と塗りつぶしの色を変更
HG4= HG3 + scale_color_manual(values = c("red", "black")) + scale_fill_manual(values = c("red", "grey")) plot(HG4)
背景の操作1
HG5= HG4 + theme_bw() plot(HG5)
背景の操作2
HG6= HG4 + theme_void() plot(HG6)
背景の操作3
HG7= HG4 + theme_minimal() plot(HG7)
背景の操作4
HG8= HG4 + theme_classic() plot(HG8)
t分布からの信頼区間
投稿: 2017-06-12
更新: 2022-1-16
例）p=0.025、自由度5の場合
# 分位点関数 quantile (q <- qt(0.025, 5, lower.tail = F)) # 確率関数 probability (p <- pt(q, 5, lower.tail = F)) # 密度関数 density (d <- dt(q, 5))
サンプルサイズ６の95％信頼区間を求める
#データサンプル DBP<-c(26,19,28,24,29,32) vDBP<-var(DBP)
統計量T値の求め方
(t <- mean(DBP) /sqrt(vDBP/6))
95%信頼区間の求め方
mean(DBP)-qt(0.025,5,lower.tail = F)*sqrt(vDBP/6) mean(DBP)+qt(0.025,5,lower.tail = F)*sqrt(vDBP/6)
Rの関数で確認
t.test(DBP) One Sample t-test data: DBP t = 14.328, df = 5, p-value = 2.985e-05 alternative hypothesis: true mean is not equal to 0 95 percent confidence interval: 21.60893 31.05774 sample estimates: mean of x 26.33333
t分布のグラフ
x <- seq(-4, 4, 0.1) plot(x, dt(x, 20), type="l") for(i in c(2,3,4)){ for(i2 in 2:4){ curve(dt(x, 6-i), col=i, type="l",add=T, ) } } labels<-c("自由度","黒:20","赤:4","緑:3","青:2") legend("topleft", legend = labels)
