中心極限定理

投稿日: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))

参考) 西内 啓; 統計学が最強の学問である[実践編]---データ分析のための思想と方法, ダイヤモンド社, 2014

Nagelkerke's R squared

Nagelkerkes's R² の求め方

logL: 最大対数尤度
deviance ＝-2*logL

Dnull: Null modelのdeviance
D: 予測したモデル（変数あり）のdeviance
n: サンプルサイズ
R² : Nagelkerkes's R²

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

参考

yoshida931.hatenablog.com

次に最尤推定法による単回帰分析（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)

CSVファイルにして作業ディレクトリに保存

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)

対応のあるデータのグラフ

id <- 1:5
x <- c(rnorm(5, 20, 2))
y <- c(rnorm(5, 25, 1))
data <- data.frame(id, x, y)

  id        x        y
1  1 22.83957 23.13902
2  2 15.05851 23.15993
3  3 19.49754 24.70116
4  4 17.43654 26.60801
5  5 17.01846 25.09940

x <- c(1,2)
t_data <- t(data[,2:3])

      [,1]     [,2]     [,3]     [,4]     [,5]
x 22.83957 15.05851 19.49754 17.43654 17.01846
y 23.13902 23.15993 24.70116 26.60801 25.09940

matplot(x, t_data, type = "l", lty = 1)
legend("topleft", legend=id, col = id,  lty = 1)

id <- 1:20
group <- c(rep("A",10),rep("B",10))
x <- c(rnorm(10,20,2), rnorm(10,15,3))
y <- c(rnorm(10,25,1), rnorm(10,20,0.5))
data <- data.frame(id, group, x, y)

x <- c(1,2)
t_data <- t(data[,3:4])
mg <- ifelse(data$group == "A", 1, 2)
matplot(x, t_data, type = "l", lty = 1, pch = 1, col = mg, xlab = "", ylab = "",xaxt="n", xlim = c(0.9,2.1))
name<-c("前","後")
axis(side=1,at=c(1, 2),labels=name)　 #指定した場所に挿入
legend(1.8, 15, c("治療A","治療B"), col = c(1, 2), lty = 1)