1級対応_統計学7章

1 問7.1.

1.1 (1) 等分散を仮定した２標本t検定.

等分散が満たされない時はwelchになる。

模擬データ.

library(MASS)
mu <- c(33.3,30.4)#平均
R <- matrix(c(1,0.7,#相関行列
              0.7,1), 2, 2)
A <- matrix(c(sqrt(36),0,
              0,sqrt(40)), 2, 2)
Sigma<- A %*% R %*% A #分散共分散行列
R#相関行列

##      [,1] [,2]
## [1,]  1.0  0.7
## [2,]  0.7  1.0

Sigma#分散共分散行列

##          [,1]     [,2]
## [1,] 36.00000 26.56313
## [2,] 26.56313 40.00000

data <- mvrnorm(20, mu, Sigma)#20個の乱数を生成
d <- data.frame(処置前=data[,1],処置後=data[,2])#データフレーム化
round(cov(d), digits = 2)#分散共分散行列

##        処置前 処置後
## 処置前  31.73  26.67
## 処置後  26.67  39.35

round(cor(d), digits = 2)#相関行列

##        処置前 処置後
## 処置前   1.00   0.75
## 処置後   0.75   1.00

par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(d)

boxplot(d)

dat <- cbind(id=rep(c(1:20),1),d)
str(dat)

## 'data.frame':    20 obs. of  3 variables:
##  $ id    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ 処置前: num  23.4 35.2 32.2 29.6 31.3 ...
##  $ 処置後: num  24.7 31.3 29.4 22.5 26 ...

1.1.1 2標本t検定

gather関数とspread関数を使いこなす.

d.gathered <- tidyr::gather(data = dat, key = 処置, value = value,処置前,処置後)
d.gathered

id	処置	value
1	処置前	23.40763
2	処置前	35.15723
3	処置前	32.17413
4	処置前	29.64306
5	処置前	31.26477
6	処置前	31.29746
7	処置前	44.96947
8	処置前	33.35679
9	処置前	36.43517
10	処置前	31.97282
11	処置前	38.25807
12	処置前	32.56909
13	処置前	28.91439
14	処置前	38.91725
15	処置前	30.21265
16	処置前	29.39011
17	処置前	41.83628
18	処置前	25.06079
19	処置前	27.71804
20	処置前	40.95700
1	処置後	24.70632
2	処置後	31.28387
3	処置後	29.36185
4	処置後	22.46479
5	処置後	26.04580
6	処置後	24.48399
7	処置後	39.29569
8	処置後	37.06753
9	処置後	35.16696
10	処置後	26.75422
11	処置後	31.05255
12	処置後	35.55926
13	処置後	34.92977
14	処置後	30.98872
15	処置後	26.45822
16	処置後	24.76077
17	処置後	39.82637
18	処置後	24.67714
19	処置後	27.69289
20	処置後	45.14646

#library("car")
t.test(d.gathered$value ~ d.gathered$処置,var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  d.gathered$value by d.gathered$処置
## t = 1.2144, df = 38, p-value = 0.2321
## alternative hypothesis: true difference in means between group 処置前 and group 処置後 is not equal to 0
## 95 percent confidence interval:
##  -1.526987  6.105890
## sample estimates:
## mean in group 処置前 mean in group 処置後 
##             33.17561             30.88616

1.1.2 手計算.

処置前と処置後が無相関を前提としている。

#プールした分散
v <- (19*var(d[,1])+19*var(d[,2]))/38
v

## [1] 35.54074

T <- mean(d[,1])-mean(d[,2])
t <- T/(sqrt(v/20*2));t

## [1] 1.214418

1.2 (2) １標本t検定.

x <- d[,1]-d[,2]
t <- (mean(x)-0)/sqrt(var(x)/20);t

## [1] 2.431377

testResult <- t.test(x, mu = 0)
testResult

## 
##  One Sample t-test
## 
## data:  x
## t = 2.4314, df = 19, p-value = 0.02511
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3186028 4.2602998
## sample estimates:
## mean of x 
##  2.289451

res<- t.test(d[,1],d[,2],paired=TRUE,alternative ="two.sided")
res

## 
##  Paired t-test
## 
## data:  d[, 1] and d[, 2]
## t = 2.4314, df = 19, p-value = 0.02511
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.3186028 4.2602998
## sample estimates:
## mean difference 
##        2.289451

1.2.1 手計算.

v <- var(d[,1])+var(d[,2])-2*cov(d[,1],d[,2])
v

## [1] 17.73324

T <- mean(d[,1])-mean(d[,2])
t <- T/(sqrt(v/20));t

## [1] 2.431377

1.2.2 可視化.

par(family = "HiraKakuProN-W3") #日本語フォントの指定
da <- cbind(id=rep(c(1:20),1),d)
da

id	処置前	処置後
1	23.40763	24.70632
2	35.15723	31.28387
3	32.17413	29.36185
4	29.64306	22.46479
5	31.26477	26.04580
6	31.29746	24.48399
7	44.96947	39.29569
8	33.35679	37.06753
9	36.43517	35.16696
10	31.97282	26.75422
11	38.25807	31.05255
12	32.56909	35.55926
13	28.91439	34.92977
14	38.91725	30.98872
15	30.21265	26.45822
16	29.39011	24.76077
17	41.83628	39.82637
18	25.06079	24.67714
19	27.71804	27.69289
20	40.95700	45.14646

#表の準備
t(da[,2:3])#転置

##            [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
## 処置前 23.40763 35.15723 32.17413 29.64306 31.26477 31.29746 44.96947 33.35679
## 処置後 24.70632 31.28387 29.36185 22.46479 26.04580 24.48399 39.29569 37.06753
##            [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]    [,16]
## 処置前 36.43517 31.97282 38.25807 32.56909 28.91439 38.91725 30.21265 29.39011
## 処置後 35.16696 26.75422 31.05255 35.55926 34.92977 30.98872 26.45822 24.76077
##           [,17]    [,18]    [,19]    [,20]
## 処置前 41.83628 25.06079 27.71804 40.95700
## 処置後 39.82637 24.67714 27.69289 45.14646

dat2 <- t(da[,2:3])
x <- c(1.1, 1.9)
matplot(
    x, dat2,
    type ="b", lty=2,
    xaxt="n", 
    xlim=c(1, 2),
    xlab="", ylab=""
)
name <- c("処置前", "処置後")
axis(side=1, at=c(1.1, 1.9), labels=name)

参考サイト.

https://sites.google.com/view/s-inf-datasci/rでt検定.
https://y2pt.com/5233/.
https://biolab.sakura.ne.jp/paired-test.html.
https://www2.kpu.ac.jp/for_ecol/obenkyou/GLMMexample.pdf.
混合モデル.

https://qiita.com/yasainiki/items/3bb65a98a8d9b0efbe64

https://rstudio-pubs-static.s3.amazonaws.com/830500_9ac256db2649449cb813798837d0604a.html#rにおけるデータの結合

2 問7.2.

2.1 上位datA

library(MASS)
mu <- c(55.89,54.74)#平均
R <- matrix(c(1,0.278,#相関行列
              0.278,1), 2, 2)
A <- matrix(c(7.21,0,
              0,12.33), 2, 2)
Sigma<- A %*% R %*% A #分散共分散行列
data <- mvrnorm(53, mu, Sigma)#20個の乱数を生成
head(data)

##          [,1]     [,2]
## [1,] 56.36112 60.40738
## [2,] 54.02866 59.90619
## [3,] 49.34616 66.81659
## [4,] 49.43449 41.79233
## [5,] 59.09370 62.43249
## [6,] 64.53334 72.25090

cov(data)

##          [,1]      [,2]
## [1,] 42.35864  20.01726
## [2,] 20.01726 106.32565

class(data)

## [1] "matrix" "array"

datA <- data.frame(class=c(rep("A",53)),pre=data[,1], post=data[,2])
datA$class <- as.factor(datA$class)
head(datA)

class	pre	post
A	56.36112	60.40738
A	54.02866	59.90619
A	49.34616	66.81659
A	49.43449	41.79233
A	59.09370	62.43249
A	64.53334	72.25090

modelA <- lm(post ~ pre ,data=datA)
summary(modelA)

## 
## Call:
## lm(formula = post ~ pre, data = datA)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -22.1522  -7.3630   0.0899   6.5580  21.7246 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  29.7778    11.5502   2.578   0.0129 *
## pre           0.4726     0.2118   2.232   0.0301 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.938 on 51 degrees of freedom
## Multiple R-squared:  0.08897,    Adjusted R-squared:  0.0711 
## F-statistic:  4.98 on 1 and 51 DF,  p-value: 0.03006

2.2 下位datB

library(MASS)
mu <- c(33.6,37.62)#平均
R <- matrix(c(1,0.304,#相関行列
              0.304,1), 2, 2)
A <- matrix(c(8.01,0,
              0,14.3), 2, 2)
Sigma<- A %*% R %*% A #分散共分散行列
data <- mvrnorm(47, mu, Sigma)#20個の乱数を生成
head(data)

##          [,1]     [,2]
## [1,] 45.98480 60.78359
## [2,] 29.96057 26.79958
## [3,] 29.98287 47.74109
## [4,] 37.48254 54.84149
## [5,] 19.87236 32.46332
## [6,] 36.51598 48.24276

cov(data)

##          [,1]      [,2]
## [1,] 80.41636  50.18012
## [2,] 50.18012 171.35775

class(data)

## [1] "matrix" "array"

datB <- data.frame(class=c(rep("B",47)),pre=data[,1], post=data[,2])
datB$class <- as.factor(datB$class)
head(datB)

class	pre	post
B	45.98480	60.78359
B	29.96057	26.79958
B	29.98287	47.74109
B	37.48254	54.84149
B	19.87236	32.46332
B	36.51598	48.24276

modelB <- lm(post ~ pre ,data=datB)
summary(modelB)

## 
## Call:
## lm(formula = post ~ pre, data = datB)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.075  -8.741  -2.221  10.322  21.642 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  16.3373     6.6088   2.472  0.01728 * 
## pre           0.6240     0.1967   3.172  0.00273 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.96 on 45 degrees of freedom
## Multiple R-squared:  0.1827, Adjusted R-squared:  0.1646 
## F-statistic: 10.06 on 1 and 45 DF,  p-value: 0.002727

2.3 データを結合

str(datA)

## 'data.frame':    53 obs. of  3 variables:
##  $ class: Factor w/ 1 level "A": 1 1 1 1 1 1 1 1 1 1 ...
##  $ pre  : num  56.4 54 49.3 49.4 59.1 ...
##  $ post : num  60.4 59.9 66.8 41.8 62.4 ...

str(datB)

## 'data.frame':    47 obs. of  3 variables:
##  $ class: Factor w/ 1 level "B": 1 1 1 1 1 1 1 1 1 1 ...
##  $ pre  : num  46 30 30 37.5 19.9 ...
##  $ post : num  60.8 26.8 47.7 54.8 32.5 ...

dat <-rbind(datA,datB)
str(dat)

## 'data.frame':    100 obs. of  3 variables:
##  $ class: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
##  $ pre  : num  56.4 54 49.3 49.4 59.1 ...
##  $ post : num  60.4 59.9 66.8 41.8 62.4 ...

2.3.1 散布図

par(family = "HiraKakuProN-W3") #日本語フォントの指定
fig1 <- function() 
    {
    pchAB <- ifelse(
        dat$class== "A", 19, 21
        )
    plot(
        dat$pre, dat$post,
        pch=pchAB, cex=1.5 ,
        xlab="入学時",  ylab="学期末"
        )
    legend(
        "topleft",
        legend = c("A", "B"),
        pch = c(19, 21)
        )
    }
fig1()

2.4 A,Bの区別をしないとき

fit <- lm(post ~ pre ,  data=dat)
summary(fit)

## 
## Call:
## lm(formula = post ~ pre, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -28.8004  -8.3255  -0.0341   8.3672  22.9588 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.89255    3.79576   3.397 0.000987 ***
## pre          0.76560    0.08269   9.259 4.91e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11 on 98 degrees of freedom
## Multiple R-squared:  0.4666, Adjusted R-squared:  0.4612 
## F-statistic: 85.73 on 1 and 98 DF,  p-value: 4.911e-15

confint(fit)

##                 2.5 %     97.5 %
## (Intercept) 5.3599950 20.4251093
## pre         0.6015079  0.9296846

attach(dat)
plot(pre,post)
abline(fit)

2.5 分散分析

#作図
par(family = "HiraKakuProN-W3") #日本語フォントの指定
fig3 <- function() 
    {
    fig1()
    
    lines(
        range(datA$pre),
        modelA$coef[1]+modelA$coef[2]*range(datA$pre),
        col="red"
     )
 
    lines(
        range(datB$pre),
        modelB$coef[1]+modelB$coef[2]*range(datB$pre),
        col="red"
    )
}
 
fig3()

2.6 共分散分析

fit2 <- lm(post ~ class*pre,  data=dat)
summary(fit2)

## 
## Call:
## lm(formula = post ~ class * pre, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.0748  -7.8548  -0.1094   8.7197  21.7246 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  29.7778    12.7089   2.343   0.0212 *
## classB      -13.4405    14.0711  -0.955   0.3419  
## pre           0.4726     0.2330   2.028   0.0453 *
## classB:pre    0.1514     0.2943   0.515   0.6080  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.93 on 96 degrees of freedom
## Multiple R-squared:  0.4837, Adjusted R-squared:  0.4675 
## F-statistic: 29.98 on 3 and 96 DF,  p-value: 9.207e-14

anova(fit2)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
class	1	8820.48757	8820.48757	73.7658386	0.0000000
pre	1	1900.61065	1900.61065	15.8948287	0.0001306
class:pre	1	31.66138	31.66138	0.2647845	0.6080347
Residuals	96	11479.11855	119.57415	NA	NA

http://nlp.dse.ibaraki.ac.jp/~shinnou/siryou/toukei-kentei/4-stat-test.pdf https://www.math.is.tohoku.ac.jp/~obata/student/subject/file/2022Stat/ELecture11.pdf

https://oshiete.goo.ne.jp/qa/9377186.html

3 問7.4.

https://www.yasuhisay.info/entry/20091018/1255826530.

mixture_gaussian <- function(x) {
  pi_0 <- 0.25
  ifelse(runif(1) < pi_0, rnorm(1, 70, 5), rnorm(1, 50, 5))
}

N <- 100
x <- sapply(1:N, mixture_gaussian)
plot(density(x)) 

x=c(0:100);x

##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100

p1 <- dnorm(x, mean=70, sd=10)*0.25
p2 <- dnorm(x, mean=50, sd=10)*0.75
p <- p1+p2
p1/p

##   [1] 2.048067e-06 2.501513e-06 3.055353e-06 3.731814e-06 4.558045e-06
##   [6] 5.567203e-06 6.799788e-06 8.305268e-06 1.014406e-05 1.238995e-05
##  [11] 1.513308e-05 1.848352e-05 2.257574e-05 2.757393e-05 3.367867e-05
##  [16] 4.113491e-05 5.024183e-05 6.136483e-05 7.495016e-05 9.154281e-05
##  [21] 1.118084e-04 1.365597e-04 1.667893e-04 2.037094e-04 2.488000e-04
##  [26] 3.038683e-04 3.711206e-04 4.532505e-04 5.535458e-04 6.760196e-04
##  [31] 8.255686e-04 1.008167e-03 1.231104e-03 1.503264e-03 1.835480e-03
##  [36] 2.240949e-03 2.735744e-03 3.339423e-03 4.075767e-03 4.973665e-03
##  [41] 6.068166e-03 7.401730e-03 9.025702e-03 1.100203e-02 1.340526e-02
##  [46] 1.632477e-02 1.986731e-02 2.415972e-02 2.935174e-02 3.561883e-02
##  [51] 4.316453e-02 5.222221e-02 6.305529e-02 7.595552e-02 9.123796e-02
##  [56] 1.092318e-01 1.302656e-01 1.546466e-01 1.826326e-01 2.143987e-01
##  [61] 2.500000e-01 2.893358e-01 3.321200e-01 3.778668e-01 4.258968e-01
##  [66] 4.753669e-01 5.253252e-01 5.747817e-01 6.227854e-01 6.684954e-01
##  [71] 7.112346e-01 7.505200e-01 7.860684e-01 8.177814e-01 8.457159e-01
##  [76] 8.700485e-01 8.910380e-01 9.089919e-01 9.242390e-01 9.371085e-01
##  [81] 9.479150e-01 9.569500e-01 9.644764e-01 9.707272e-01 9.759056e-01
##  [86] 9.801867e-01 9.837197e-01 9.866314e-01 9.890281e-01 9.909991e-01
##  [91] 9.926186e-01 9.939486e-01 9.950401e-01 9.959355e-01 9.966698e-01
##  [96] 9.972718e-01 9.977652e-01 9.981696e-01 9.985009e-01 9.987723e-01
## [101] 9.989946e-01

plot(x,p,type="l")

x=c(0:100);x

##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100

p1 <- dnorm(x, mean=70, sd=10);p1

##   [1] 9.134720e-13 1.830332e-12 3.630962e-12 7.131328e-12 1.386680e-11
##   [6] 2.669557e-11 5.088140e-11 9.601433e-11 1.793784e-10 3.317884e-10
##  [11] 6.075883e-10 1.101576e-09 1.977320e-09 3.513955e-09 6.182621e-09
##  [16] 1.076976e-08 1.857362e-08 3.171349e-08 5.361035e-08 8.972435e-08
##  [21] 1.486720e-07 2.438961e-07 3.961299e-07 6.369825e-07 1.014085e-06
##  [26] 1.598374e-06 2.494247e-06 3.853520e-06 5.894307e-06 8.926166e-06
##  [31] 1.338302e-05 1.986555e-05 2.919469e-05 4.247803e-05 6.119019e-05
##  [36] 8.726827e-05 1.232219e-04 1.722569e-04 2.384088e-04 3.266819e-04
##  [41] 4.431848e-04 5.952532e-04 7.915452e-04 1.042093e-03 1.358297e-03
##  [46] 1.752830e-03 2.239453e-03 2.832704e-03 3.547459e-03 4.398360e-03
##  [51] 5.399097e-03 6.561581e-03 7.895016e-03 9.404908e-03 1.109208e-02
##  [56] 1.295176e-02 1.497275e-02 1.713686e-02 1.941861e-02 2.178522e-02
##  [61] 2.419707e-02 2.660852e-02 2.896916e-02 3.122539e-02 3.332246e-02
##  [66] 3.520653e-02 3.682701e-02 3.813878e-02 3.910427e-02 3.969525e-02
##  [71] 3.989423e-02 3.969525e-02 3.910427e-02 3.813878e-02 3.682701e-02
##  [76] 3.520653e-02 3.332246e-02 3.122539e-02 2.896916e-02 2.660852e-02
##  [81] 2.419707e-02 2.178522e-02 1.941861e-02 1.713686e-02 1.497275e-02
##  [86] 1.295176e-02 1.109208e-02 9.404908e-03 7.895016e-03 6.561581e-03
##  [91] 5.399097e-03 4.398360e-03 3.547459e-03 2.832704e-03 2.239453e-03
##  [96] 1.752830e-03 1.358297e-03 1.042093e-03 7.915452e-04 5.952532e-04
## [101] 4.431848e-04

plot(x,p1)

3.1 (1)

x <- 60
p1 <- dnorm(x, mean=70, sd=10)*0.25
p2 <- dnorm(x, mean=50, sd=10)*0.75
p <- p1+p2
p1/p

## [1] 0.25

3.2 (2)

x <- 70
p1 <- dnorm(x, mean=70, sd=10)*0.25
p2 <- dnorm(x, mean=50, sd=10)*0.75
p <- p1+p2
p1/p

## [1] 0.7112346

3.3 (3)

x=c(1:100);x

##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100

p1 <- dnorm(x, mean=70, sd=10)*0.25
p2 <- dnorm(x, mean=50, sd=10)*0.75
pd <- data.frame(p=p1/(p1+p2))

subset(pd,p>0.8)

	p
73	0.8177814
74	0.8457159
75	0.8700485
76	0.8910380
77	0.9089919
78	0.9242390
79	0.9371085
80	0.9479150
81	0.9569500
82	0.9644764
83	0.9707272
84	0.9759056
85	0.9801867
86	0.9837197
87	0.9866314
88	0.9890281
89	0.9909991
90	0.9926186
91	0.9939486
92	0.9950401
93	0.9959355
94	0.9966698
95	0.9972718
96	0.9977652
97	0.9981696
98	0.9985009
99	0.9987723
100	0.9989946