1級対応_統計学_6章
1 問１
1.1 データ
ID <- c(1:10);ID## [1] 1 2 3 4 5 6 7 8 9 10L <- c(250,300,300,200,200,200,250,350,200,250);L## [1] 250 300 300 200 200 200 250 350 200 250R <- c(150,250,300,150,200,200,200,300,200,150);R## [1] 150 250 300 150 200 200 200 300 200 150Y <- c(53,74,99,66,45,69,59,86,38,46);Y## [1] 53 74 99 66 45 69 59 86 38 46data <- data.frame(ID=factor(ID),L=L,R=R,T=L+R,Y=Y)
data| ID | L | R | T | Y |
|---|---|---|---|---|
| 1 | 250 | 150 | 400 | 53 |
| 2 | 300 | 250 | 550 | 74 |
| 3 | 300 | 300 | 600 | 99 |
| 4 | 200 | 150 | 350 | 66 |
| 5 | 200 | 200 | 400 | 45 |
| 6 | 200 | 200 | 400 | 69 |
| 7 | 250 | 200 | 450 | 59 |
| 8 | 350 | 300 | 650 | 86 |
| 9 | 200 | 200 | 400 | 38 |
| 10 | 250 | 150 | 400 | 46 |
library(psych)
#describe(data[,2:5])1.2 基本事項
- 重相関係数は予測値と実測値との相関係数.
- 決定係数 は回帰変動／全変動.
R^2=予測値の分散／実測値の分散.
- F値=(回帰分散)／（残差分散）.
res <- lm(Y~L+R)
#str(res)
summary(res)##
## Call:
## lm(formula = Y ~ L + R)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.615 -5.612 -2.381 9.571 19.077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.93846 21.50225 -0.090 0.931
## L 0.09892 0.12578 0.786 0.457
## R 0.19385 0.11678 1.660 0.141
##
## Residual standard error: 13.32 on 7 degrees of freedom
## Multiple R-squared: 0.6287, Adjusted R-squared: 0.5226
## F-statistic: 5.926 on 2 and 7 DF, p-value: 0.03119fit <- res $ fitted.values
#重相関係数は予測値と実測値との相関係数
cor(fit,Y)## [1] 0.792899#決定係数は回帰変動／全変動
#R^2=予測値の分散／実測値の分散=(SR/9)/(ST/9)
var(fit)/var(Y)## [1] 0.6286888#回帰変動
var(fit)*9## [1] 2101.392#全変動
var(Y)*9## [1] 3342.5(9*var(Y)-9*var(fit))/7#残差分散## [1] 177.3011#F値
FF <- var(fit)*9/2/((9*var(Y)-9*var(fit))/7);FF## [1] 5.926056#F値=(回帰分散)/（残差分散）
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:psych':
##
## logit#各変数に多重共線性が無いかを調べるための統計量。 各変数の相関係数も多重共線性の指標となりうる
vif( res )## L R
## 2.230769 2.230769- 実測値の変動から予測値の変動を引いたものは誤差変動である。
それを誤差の自由度で割れば誤差分散になる。 - 予測値の変動を変数の数で割った回帰分散は、回帰がないという帰無仮説のもとでは誤差分散と一致する。即ち、誤差の分しか変動しない。しかし、回帰があれば回帰分散は大きくなりF値も大きくなる。
fit = lm( Y~L+R, data = data )
library(car)
vif( fit )## L R
## 2.230769 2.230769fit = lm( Y~L+R+L:R, data = data )
summary( fit )##
## Call:
## lm(formula = Y ~ L + R + L:R, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.209 -9.161 1.603 9.428 14.490
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 126.428082 157.831714 0.801 0.454
## L -0.390479 0.609618 -0.641 0.545
## R -0.383356 0.712856 -0.538 0.610
## L:R 0.002114 0.002573 0.821 0.443
##
## Residual standard error: 13.64 on 6 degrees of freedom
## Multiple R-squared: 0.6662, Adjusted R-squared: 0.4993
## F-statistic: 3.992 on 3 and 6 DF, p-value: 0.07037library( effects )## lattice theme set by effectsTheme()
## See ?effectsTheme for details.plot( effect( "L:R", fit ), multiline = T)
多重共線性は VIF (variance inflation factor; 分散拡大要因) という指標で診断できます。 VIF は各説明変数に対して算出され、その値は 2 未満であることが望まれます。 VIF の値が 10 以上であれば、その変数は多重共線性の要因になっていることが示唆されます。
「2_回帰の多重共線性」 参照
2 問２
library(MASS)
mu <- c(58.4,48.917,44.033,49.283,59.508)#平均
R <- matrix(c(1,0.84,0.206,0.283,0.577,#相関行列
0.84,1,0.366,0.426,0.62,
0.206,0.366,1,0.819,0.453,
0.283,0.426,0.819,1,0.464,
0.577,0.62,0.453,0.464,1), 5, 5)
A <- matrix(c(13.81,0,0,0,0,
0,9.772,0,0,0,
0,0,13.769,0,0,
0,0,0,7.514,0,
0,0,0,0,7.503), 5, 5)
Sigma<- A %*% R %*% A #分散共分散行列
data <- mvrnorm(120, mu, Sigma)#120個の乱数を生成
d <- data.frame(国=data[,1],社=data[,2],数=data[,3],理=data[,4],英=data[,5])#データフレーム化
round(cov(d), digits = 2)#分散共分散行列## 国 社 数 理 英
## 国 199.55 129.34 44.95 34.63 70.31
## 社 129.34 111.05 48.88 31.20 55.47
## 数 44.95 48.88 164.44 73.22 49.60
## 理 34.63 31.20 73.22 48.69 29.37
## 英 70.31 55.47 49.60 29.37 63.50round(cor(d), digits = 2)#相関行列## 国 社 数 理 英
## 国 1.00 0.87 0.25 0.35 0.62
## 社 0.87 1.00 0.36 0.42 0.66
## 数 0.25 0.36 1.00 0.82 0.49
## 理 0.35 0.42 0.82 1.00 0.53
## 英 0.62 0.66 0.49 0.53 1.002.1 標準化したデータ
data <- mvrnorm(120, mu, R)#120個の乱数を生成
d <- data.frame(国=data[,1],社=data[,2],数=data[,3],理=data[,4],英=data[,5])#データフレーム化
round(cov(d), digits = 2)#分散共分散行列## 国 社 数 理 英
## 国 1.29 1.18 0.23 0.41 0.65
## 社 1.18 1.40 0.40 0.62 0.68
## 数 0.23 0.40 1.12 0.95 0.53
## 理 0.41 0.62 0.95 1.17 0.63
## 英 0.65 0.68 0.53 0.63 1.02round(cor(d), digits = 2)#相関行列## 国 社 数 理 英
## 国 1.00 0.88 0.19 0.34 0.57
## 社 0.88 1.00 0.32 0.48 0.57
## 数 0.19 0.32 1.00 0.82 0.50
## 理 0.34 0.48 0.82 1.00 0.57
## 英 0.57 0.57 0.50 0.57 1.002.2 主成分分析（分散行列）
library(psych)
summary(d)## 国 社 数 理
## Min. :55.22 Min. :45.37 Min. :41.37 Min. :46.06
## 1st Qu.:57.56 1st Qu.:48.06 1st Qu.:43.26 1st Qu.:48.45
## Median :58.39 Median :48.90 Median :43.99 Median :49.07
## Mean :58.39 Mean :48.80 Mean :43.95 Mean :49.14
## 3rd Qu.:59.26 3rd Qu.:49.57 3rd Qu.:44.66 3rd Qu.:49.78
## Max. :61.08 Max. :51.71 Max. :47.02 Max. :51.73
## 英
## Min. :57.30
## 1st Qu.:58.81
## Median :59.37
## Mean :59.42
## 3rd Qu.:60.08
## Max. :61.89result <- prcomp(d, scale=F) #scale=Tは相関行列から分析
summary(result)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.938 1.213 0.67558 0.43279 0.37115
## Proportion of Variance 0.625 0.245 0.07594 0.03116 0.02292
## Cumulative Proportion 0.625 0.870 0.94592 0.97708 1.000002.3 主成分分析（相関行列）
result <- prcomp(d, scale=T)
eigen_r <- eigen(cor(d))#相関行列の固有値
round(eigen_r$vectors,digits=2)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.43 0.54 0.15 0.23 0.66
## [2,] 0.47 0.41 0.34 -0.08 -0.69
## [3,] 0.40 -0.58 0.19 0.68 -0.06
## [4,] 0.46 -0.44 0.23 -0.69 0.25
## [5,] 0.46 0.01 -0.88 -0.03 -0.11par(family = "HiraKakuProN-W3") #日本語フォントの指定
biplot(result)
2.4 負荷量
library(psych)
dat <- d
fit <- principal(dat, nfactors=2, rotate=FALSE)#回転ない## Specified rotation not found, rotate='none' used#str(fit)
fit $loadings##
## Loadings:
## PC1 PC2
## 国 0.761 -0.598
## 社 0.832 -0.454
## 数 0.708 0.640
## 理 0.813 0.482
## 英 0.818
##
## PC1 PC2
## SS loadings 3.103 1.205
## Proportion Var 0.621 0.241
## Cumulative Var 0.621 0.862par(family = "HiraKakuProN-W3") #日本語フォントの指定
biplot(fit)
3 問３
文12_共1 の類題 「回帰分析
決定係数、予測値の信頼区間」
高校ユニット「回帰分析」参照
X <- c(0,1,2,3,4)
Y1 <- c(0.5,0,1.5,1,2.5)
data1 <- data.frame(X=X,Y1=Y1);data1| X | Y1 |
|---|---|
| 0 | 0.5 |
| 1 | 0.0 |
| 2 | 1.5 |
| 3 | 1.0 |
| 4 | 2.5 |
plot(data1,ylim=c(0,3))
fit1 = lm( formula = Y1 ~ X, data = data1 )
fity1 <- fit1 $ fitted.values;fity1## 1 2 3 4 5
## 0.1 0.6 1.1 1.6 2.1plot(data1,ylim=c(0,3))
abline(fit1)
a <- fit1$ coefficients;a#回帰係数## (Intercept) X
## 0.1 0.5e <- fit1$ residuals#残差
se <- sum(e^2);se#残差平方和## [1] 1.2#se/2#残差分散
sqrt(se/3)#残差の標準誤差## [1] 0.6324555#回帰係数を手計算で求める
cov(X,Y1)/var(X)## [1] 0.5#決定係数
var(fity1)/var(Y1)## [1] 0.67567574 問４
par(mfrow = c(1,2)) # グラフの描画範囲を四分割
n <- 96
d <- c(rep(0,35),rep(1,35),rep(2,18),rep(3,5),rep(4,2),rep(7,1))
lambda <- mean(d);lambda## [1] 1.052083#理論分布
x <- rpois(96,lambda)
hist(x,prob=T,breaks = c(0:5))
plot(0:5,dpois(0:5,lambda))
x## [1] 1 0 1 4 1 2 2 4 0 1 0 2 1 1 1 1 1 1 2 3 1 2 0 0 1 1 0 1 1 2 1 1 1 1 2 1 3 2
## [39] 0 3 0 0 1 1 1 2 2 1 2 0 1 0 1 1 2 1 0 1 1 2 1 1 0 0 1 0 0 2 0 0 1 1 1 3 1 1
## [77] 0 1 0 0 0 2 0 2 1 0 0 2 1 1 3 2 1 2 0 44.1 適合度検定
dpois(0:4,lambda)## [1] 0.34920947 0.36739746 0.19326637 0.06777744 0.0178268896*dpois(0:4,lambda)## [1] 33.524109 35.270157 18.553572 6.506635 1.711380#期待度数が５以下のセルをまとめる
p4 <- 1-sum(dpois(0:3,lambda))
p=c(dpois(0:3,lambda),p4)
npi <- 96*p
ni <- c(35,35,18,5,3)
(v <- (ni - npi)^2 / npi)## [1] 0.064975734 0.002069301 0.016516599 0.348866653 0.340299966# 和をとって、自由度3のカイ二乗分布の累積確率をもとに検定
(ch2 <- sum(v))## [1] 0.77272831 - pchisq(ch2, 5 - 1 - 1)#自由度に注意## [1] 0.8559783(chisq.test(ni,p=p))## Warning in chisq.test(ni, p = p): Chi-squared approximation may be incorrect##
## Chi-squared test for given probabilities
##
## data: ni
## X-squared = 0.77273, df = 4, p-value = 0.9421#カイ二乗値は同じだが自由度が異なるのでP値が異なる
vi <-chisq.test(ni,p=p)$ statistic## Warning in chisq.test(ni, p = p): Chi-squared approximation may be incorrect1-pchisq(vi, df=3)## X-squared
## 0.8559783