いろいろな統計分析手法 〜μ'sを例にして〜

統計分析の授業で習った手法を使ってμ'sを分析していきます。Rの使い方をまとめた自分用のメモなので解析手法が妥当かどうかなどは一切考えていません。

データの準備と読み込み

ラブライブ!というアニメを観ていてにあったμ'sのデータを利用します。

	B	W	H	height	age	
東條希	90	60	82	156	17
絢瀬絵里	88	60	84	162	17
小泉花陽	82	60	83	156	15
南ことり	80	58	80	159	16
西木野真姫	78	56	83	161	15
高坂穂乃果	78	58	82	157	16
園田海未	76	58	80	159	16
星空凛	75	59	80	155	15
矢澤にこ	74	57	79	154	17

これをlovelove.txtという名前で保存します。

まずはデータを読み込みます。

data <- read.table("lovelive.txt", header=T)
data

出力

            B  W  H height age
東條希     90 60 82    156  17
絢瀬絵里   88 60 84    162  17
小泉花陽   82 60 83    156  15
南ことり   80 58 80    159  16
西木野真姫 78 56 83    161  15
高坂穂乃果 78 58 82    157  16
園田海未   76 58 80    159  16
星空凛     75 59 80    155  15
矢澤にこ   74 57 79    154  17

これでデータを読み込みます。
この書き方は苗字と名前の間にスペースが入っているとうまく動かないようです。

クラスター分析

クラスター分析は複数のデータのうち近いものはどれかを調べるような解析手法です。
ここではマンハッタン距離を使った群均等法によるクラスター分析を行っています。

data.dist <- dist(data,"manhattan")
#data.dist で距離行列の出力
data.clust <- hclust(data.dist,"average")
plot(data.clust, main="μ'sのクラスター分析", sub="", xlab="名前", ylab=" 成分量の差")

結果

なんとなく近い学年と似た体型っぽい人によるクラスターができているように見えます。

対応分析

対応分析はアンケートで似た解答パターンをしている人の分析などに使う手法で、このデータの分析には向きませんがそのままやってみます

library(MASS)
corresp(data) -> data.corresp
plot(data.corresp)

この図からはあまり違いを読み取ることができません……

corresp(data, nf=2) -> data.corresp
plot(data.corresp)

Bに特徴のある人がいるなぁくらいのことは分かります。

主成分分析

主成分分析は3次元以上のデータを2次元で可視化するときなどに使う手法で、なるべく投影したときのばらつきが大きくなるような投影が行われます。
まずは各主成分の寄与率を求めます

data.pca <- prcomp(data)
plot(data.pca)
summary(data.pca)

                          PC1    PC2     PC3     PC4     PC5
Standard deviation     5.9287 2.7019 1.20465 0.85421 0.53456
Proportion of Variance 0.7825 0.1625 0.03231 0.01625 0.00636
Cumulative Proportion  0.7825 0.9451 0.97739 0.99364 1.00000

第一主成分までで78%, 第二主成分までで94%説明できるようです。

第一・第二主成分を用いた平面グラフを書いてみます。

biplot(data.pca,choices=1:2)

第一主成分は主にスリーサイズが影響していることが分かります。第二主成分は身長の影響が若干強いことが分かりますがなんだかよく分かりません。（この分析、元ネタとほぼ同じじゃないか……）

（以下授業の進行につれて更新）

プロットの保存

PNGにするとき

png("cluster.png",height=960, width=960, res=144)
plot(......)
dev.off()

PDFにするとき

pdf(file="pca.pdf", family="Japan1")
plot(.....)
dev.off()

フォントのレンダリングが微妙だから高解像度のPNGの方がよさそう

結論

どんな統計手法もうまく解釈しないと役に立たない

小テスト問題と解答

データ exercise1.txt をダウンロードしてRで読み込み、zが“A”の行のみを抽出しなさい。こうして得られた部分データの30行目のxの値とyの値の積を求めなさい。

exercise1 <- read.table("exercise1.txt",header=T)
u <- exercise1[,"z"]=="A" # データの変量zがAか否か
v <- exercise1[u,] # u がTRUEの行を抽出し、vにしまう。
w <- v[30,] # vの30行目
prod(w[1:2])

授業で用いた大腸がんの患者さんの遺伝子発現データcolon.txtを読み込み、このデータの1番目の患者さんの癌組織における遺伝子の発現量と正常組織における発現量の比の平均を求めなさい。

colon <- read.table("colon.txt",header=T)
tn.ratio <- colon[,1]/colon[,2]
mean(tn.ratio)

授業で用いた「長寿社会を支える家族・コミュニティ・制度」調査:目黒区データ(meguro96.csv)を読み込み、問2.4 (最近、自分と親との間の考え方の違いで、もめたり悩んだりしたことがありますか?)について、女性20歳代の回答パターンを分析しなさい。2.(考え方の違いはあるが、うまくやっている)を選んだ人は何パーセントか、求めなさい。

read.csv("meguro96.csv",header=T)->meguro
v <- meguro$q5.1==2 & meguro$q5.2==1 # 性別(問5.1)が女性(2)で、かつ年齢(問5.2)が20歳代(1)か否か
parent_F20 <- meguro[v,"q2.4"] # 女性20歳代の問2.4への回答データ
table(parent_F20) # 単純集計:度数表
prop.table(table(parent_F20))*100 # 単純集計:パーセント

データ exercise2.txt をマンハッタン距離・群均等法でクラスター分析しなさい。
グループ(S1,S7)とグループ(S2,S3,S5)の間の距離はどの程度か。

read.table("exercise2.txt",header=T)->exercise
plot(hclust(dist(exercise,"manhattan"),"average"))

データ exercise2.txtに数量化III類（対応分析）を適用しなさい。次の項目の組み合わせのうち、グループ(S2,S3,S5)で大きな値をとり、このグループを特徴づけているものは
1.(A,B), 2.(A,C), 3.(A,E), 4.(A,F), 5.(A,G)　のうちどれか。ひとつ選びなさい。

library(MASS)
plot(corresp(exercise))
plot(corresp(exercise,nf=2))

データrelation.txtをマンハッタン距離・群均等法でクラスター分析しなさい。男性と女性で同じようなカテゴリーを選択していることが伺われる年齢層はどれか。

read.table("relation.txt",header=T)->relation
relation.dist <- dist(relation,"manhattan")
relation.clust <- hclust(relation.dist,"average")
plot(relation.clust,main="",sub="",xlab="")

データrelation.txtに数量化III類（２次元）を適用しなさい。女性40歳代を特徴づけることが伺われるカテゴリーはどれか。

plot(corresp(relation,nf=2),xlab="1st xis",ylab="2nd axis")

データrat.txtは、30頭のラットR1‐R30について、生後8日、15日、22日、29日、36日における体重(g)を測定したものである。このデータを主成分分析しなさい。この結果から、R2, R9, R11, R14, R27のうちどの個体が最も速く成長したと推察されるか。

read.table("rat.txt",header=T)->rat
prcomp(rat)->rat.pca
biplot(rat.pca)
names(rat.pca)
rat.pca$rotation

同じくデータrat.txtの主成分分析の結果から、R1, R6, R10, R15, R25のうちどの個体がこの5週間の平均体重が最も小さかったと推察されるか。

上に同じ

大腸がんの患者さんの正常組織と癌組織における遺伝子発現のデータcolon.txtを読み込み、データを対数変換したものを主成分分析しなさい。累積寄与率が80%を超えるには、少なくともいくつの主成分が必要か。

read.table("colon.txt",header=T) -> colon
colon_log <- log(colon)
prcomp(colon_log)-> colon_log.pca
plot(colon_log.pca)
summary(colon_log.pca)

同じくcolon.txtを対数変換したデータで、第1主成分が大きい遺伝子はどのような遺伝子か。

biplot(colon_log.pca)

平均0.5, 標準偏差0.2の正規分布の上側30%点を求めなさい。

qnorm(0.3,mean=0.5,sd=0.2,lower.tail=F)

ある変数は平均0.5, 標準偏差0.2の正規分布に従うとする。この変数が0.7以上の値をとる確率を求めなさい。

pnorm(0.7,mean=0.5,sd=0.2,lower.tail=F)

添付のexam4.txtは10000個の数字からなるデータである。このデータをダウンロードして読み込みなさい。これを母集団とみなす。母平均を求めなさい

scan("exam4.txt")->pop
mean(pop)

問1の母集団から単純ランダムサンプリングにより得られた大きさ100の標本に基づき、標本平均により母平均を推定する。シミュレーションにより、この推定量の標準誤差を計算しなさい。

n.sample <- 100
trial <- 10000
samp.mean_set <- NULL
for (i in 1:trial){
samp <- sample(pop, n.sample)
samp.mean <- mean(samp)
c( samp.mean ,samp.mean_set) -> samp.mean_set
}
sd(samp.mean_set)

データexercise5.txtは、母平均を推定するために、母集団から単純ランダムサンプリングにより得られた大きさ50の標本である。このデータをダウンロードして読み込み、標本平均を計算しなさい。

scan("exercise5.txt") -> x
mean(x) -> samp.mean
samp.mean

問1で得られた標本平均の標準誤差を計算しなさい。

sd(x)/sqrt(50)

データspeak_sex.txtは近所づきあいの有無(有:1, 無: 0)と性別(男: M, 女: F)への回答を表している。このデータを読み込み、女性で近所づきあい有と答えた人の割合と男性で近所づきあい有と答えた人の割合の差を計算しなさい。

read.table("speak_sex.txt",header=T) -> speak_sex
head(speak_sex)
speak_f <- speak_sex[speak_sex$sex=="F",1]
mean(speak_f)
speak_m <- speak_sex[speak_sex$sex=="M",1]
mean(speak_m)
mean(speak_f)-mean(speak_m)

別解

speak_sex.table <- prop.table(table(speak_sex),2)
speak_sex.table
speak_sex.table[2,1] - speak_sex.table[2,2]

近所づきあいをしている人の割合に男女差がないという帰無仮説をt検定により両側検定しなさい。p値を求めなさい。

t.test(speak_f, speak_m)
t.test(speak~sex,data=speak_sex)

データexercise6.txtを読み込んで、変量yに対し、factor_Aの一元配置分散分析を行いなさい。p値はいくつか。

read.table("exercise6.txt",header=T)->exercise
boxplot(y~factor_A,exercise,col="brown",
xlab="category of factor_A",ylab="y")
aov(y~factor_A,exercise)->aov_factorA
summary(aov_factorA)

同じくデータexercise6.txtについて、変量yに対し、factor_Aとfactor_Bの交互作用も含めた二元配置分散分析を行いなさい。factor_Bの主効果のp値はいくつか。

boxplot(y~factor_A+factor_B,exercise,col="brown",
xlab="category of factor_A and factor_B",ylab="y")
aov(y~factor_A*factor_B,exercise)->aov_factorAB
summary(aov_factorAB)

表1の分散分析表において、総平方和は530とする。空欄[a]‐[h]には数値が入る。要因Bのp値(分子の自由度3,分母の自由度hのF分布で値がd以上になる確率)を求めなさい。

自由度	平方和	平均平方和	F値
要因A	2	[a]	[b]	5.6
要因B	3	42	[c]	[d]
要因A×要因B	6	[e]	[f]	[g]
誤差	[h]	240	5

b/5= 5.6
c=42/3
a/2=28
f=192/6
240/h=5
d=14/5
56+42+e+240=530
g=32/5

pf(2.8,df1=3,df2=48,lower.tail=F)

表1で、要因Aと要因Bの交互作用のp値を求めなさい。

pf(6.4,df1=6,df2=48,lower.tail=F)

添付のデータexercise7.txtを読み込んで、変量yを従属変数、変量x1を説明変数として、単回帰分析しなさい。x1の回帰係数のp値はいくつか。

read.table("exercise7.txt",header=T)->exercise
plot(y~x1,exercise)
x1.lm <- lm(y~x1,exercise)
summary(x1.lm)

同じくデータexercise7.txtについて、変量yを従属変数、変量x1とx2を説明変数として、重回帰分析しなさい。x1の回帰係数のp値はいくつか。

plot(exercise)
x12.lm <- lm(y~x1+x2,exercise)
summary(x12.lm)

添付のデータxy_thr3.txtをプロットすると、x<3ではyは変化する傾向はなく、x>3でxとともに直線的にyも増加する傾向を持つことがわかる。これにホッケー・スティックモデルをあてはめなさい。xの値が閾値3以下のとき、yの値はいくつと推定されるか。

data0$x3 <- pmax(data0$x-3,0)
data0
lm(y~x3,data0) -> xy_thr3.lm
summary(xy_thr3.lm)
data0$x2 <- data0$x^2
lm(y~x+x2,data0) -> xy_pol2.lm
summary(xy_pol2.lm)

添付のデータxy_poisson.txtは、従属変数yが説明変数xとともに平均的に増加するカウントデータである。このデータを単回帰分析しなさい。xの回帰係数はいくつか。

read.table("xy_poisson.txt",header=T)->data0
plot(y~x,data0)
lm(y~x,data0) -> xy.lm
summary(xy.lm)
abline(xy.lm,col="red")

同じくデータxy_poisson.txtを一般化線形モデルでポアソン回帰しなさい。ただし、リンク関数はデフォルトの対数変換とします。回帰係数はいくつか。

glm(y~x,poisson,data0)->xy.poisson
summary(xy.poisson)

添付のデータinhibitor.txtはある酵素に対する阻害剤の効果を調べた実験を表している。この阻害剤は酵素に競合的に働き、Vmaxを変えることなくKmを上げることがわかっている。concは基質の濃度、inv_concはその逆数、rateは反応速度、inv_rateはその逆数、
treatedは阻害剤の投与(1)・非投与(0)を表している。inv_rateを2つの説明変数inv_conc, inv_conc0に重回帰し、阻害剤の効果を推定しなさい。最大反応速度の逆数の推定値はいくつか。

read.table("inhibitor.txt",header=T)->inhibitor
inhibitor.lm <- lm(inv_rate ~ inv_conc+inv_conc0,inhibitor)
summary(inhibitor.lm)

図1はカウントデータである従属変数yを説明変数xでポアソン回帰した様子を示している。リンク関数はデフォルトの対数変換とし、一般化線形モデルglmで分析したものである。表1は回帰係数の推定値と推定精度を示している。説明変数の値が‐1のとき、従属変数yが0となる確率の推定値を求めなさい。

dpois(0,lambda=3.261966)

添付のデータneighbor1.txtは「長寿社会を支える家族・コミュニティ・制度」調査における近所づきあいの程度を尋ねた問1.1と年齢を尋ねた問5.2への回答を表している。問1.1にはカテゴリーが4つあるが、変量relationは、1または2を選択した場合は"little"、3または4を選択した場合は"intimate"と、2つのカテゴリーにまとめている。このデータをクロス集計しなさい。40歳代で"intimate"と答えた人は何人か。

read.table("neighbor1.txt",header=T)->neighbor1
table(neighbor1)->neighbor1.tab
neighbor1.tab

問1で求められたクロス集計表から、relationの2つのカテゴリーの割合(%)を年齢別に求めなさい。30歳代の回答者のうち、intimateと答えた人の割合は何パーセントか。

table(neighbor1)->neighbor1.tab
neighbor1.tab
prop.table(neighbor1.tab,1) -> neighbor1.age
round(neighbor1.age*100,3)

授業で用いた食道がんのケースコントロールスタディーのデータesoph0.txtを読み込んで、年齢群別にケースとコントロールの数を示す二重クロス集計を求めなさい。そして、カイ二乗検定により年齢と食道がんの関連の有無を調べ、p値を求めなさい。

read.table("esoph0.txt")->esoph
attach(esoph)
case.age <- tapply(ncases,agegp,sum)
control.age <- tapply(ncontrols,agegp,sum)
esoph.agegp <- data.frame(case=case.age,control=control.age)
esoph.agegp
chisq.test(esoph.agegp)-> esoph.agegp.test
esoph.agegp.test

同じくデータesoph0.txtを分析して得られた年齢群別にケースとコントロールの数を示す二重クロス集計表から、45歳群におけるcontrolに対するcaseのオッズは25歳群のオッズの何倍になっているか、オッズ比を求めなさい。

esoph.agegp
esoph.agegp.odds <- esoph.agegp[,1]/ esoph.agegp[,2]
esoph.agegp.odds
esoph.agegp.or <- esoph.agegp.odds[-1]/esoph.agegp.odds[1]
esoph.agegp.or

添付のデータxy01.txtは毒性を測る実験の模擬データである。投与量xに対してy(1が生存、0が死亡)が記録されている。yをxに対してロジスティック回帰しなさい。生存率が50%になるときの投与量はいくつと推定されるか。

read.table("xy01.txt",header=T)->data0
xy.logist <- glm(y~x,binomial,data0)
summary(xy.logist)
#log(p/1-p)=-8.34+3.78xなのでこれが0のときのx, x=2.21が答え

同じく添付のデータxy01.txtにロジスティック回帰を行った結果を下に、x=1のときの生存率を推定しなさい。
p=1/(1+exp(-8.34+3.78x)でx=1を代入すればよい

授業で用いた食道がんのケースコントロールスタディーのデータesoph0.txtを読み込んで、症例vs対照群(ncases, ncontrols)を喫煙量の因子(tobgp)にロジスティック回帰しなさい。喫煙量のカテゴリーT0 (0‐9g/day)に対するカテゴリーT10 (10‐19g/day)のオッズ比はいくつか。

read.table("esoph0.txt",header=T)->data0
esoph.tob <- glm(cbind(ncases,ncontrols)~tobgp,binomial,data0)
summary(esoph.tob)
#対数だからexpを取る

添付のmature_age_meta0.txtは10個の局所集団(local population)からの標本について年齢(age)と性状態(maturity, 0:未成熟、1成熟)を記録した模擬データである。性状態を年齢にロジスティック回帰(family=binomial)しなさい。成熟率をロジット変換したときの切片の値を求めなさい。

read.table("mature_age_meta0.txt")->mature_age_meta
age_maturity.glm <-glm(maturity~age,family=binomial,data=mature_age_meta)
summary(age_maturity.glm)

同じくmature_age_meta0.txtを分析しなさい。ただし今度は、性状態を年齢にロジスティック回帰する際に、切片の局所集団(local population: loc_pop)の間での変動をランダム効果として取り入れた一般化線形混合モデルを当てはめなさい。成熟率をロジット変換し
たときの切片が局所集団(local population: loc_pop)の間でどの程度変動するか、標準偏差を求めなさい。

library(lme4)
age_maturity.glmm <-glmer(maturity~age+(1|loc_pop),
family=binomial,data=mature_age_meta)
summary(age_maturity.glmm)

添付の模擬データxy.txtを読みこみ、yをxに単回帰しなさい。これを下に、x=2に対する予測の95%信頼区間を求めなさい。

read.table("xy.txt",header=T)->data0
xy.lm <- lm(y~x,data0)
x.data <- data.frame(x=2)
xy.pred <- predict(xy.lm,x.data,interval="confidence")
xy.pred

同じく添付の模擬データxy.txtについて、x=5に対する95%予測区間を求めなさい。

x.data2 <- data.frame(x=5)
xy.pred2 <- predict(xy.lm,x.data,interval="prediction")
xy.pred2

添付のneighbor_e.txtを読み込み、relationを他のすべての変数に対してロジスティック回帰しなさい。こうして得られた女性に対する男性の近所づきあいのオッズ比はいくつか。

read.table("neighbor_e.txt",header=T)->neighbor1
neighbor.glm <- glm(relation~.,binomial,neighbor1)
summary(neighbor.glm)
exp(-0.608)=0.544 #対数だから戻す

問1で行ったロジスティック回帰にstepを適用して、変数選択を行いなさい。選択されたモデルのAICはいくつか。

step(neighbor.glm) -> neighbor.step
summary(neighbor.step)