ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2011-08-03 ぱらぱらめくる『Statistics for Experimenters』

[][][]ぱらぱらめくる『Statistics for Experimenters』

  • 目次
    • Ch1 Catalyzing the Generation of Knowledge
    • Ch2 Basics (Probability,Parameters, and Statistics)
    • Ch3 Comparing Two Entities: Reference Distributions, Tests, and Confidence Intervals
    • Ch4 Comparing a Number of Entities, Randomized Blocks, and Latin Squares
    • Ch5 Factorial Designs at Two Levels
    • Ch6 Fractional Factorial Designs
    • Ch7 Additional Fractionals and Analysis
    • Ch8 Factorial Designs and Data Transformation
    • Ch9 Multiple Sources of Variation
    • Ch10 Least Squares and Why We Need Designed Experiments
    • Ch11 Modeling, Geometry, and Experimental Design
    • Ch12 Some Applications of Responce Surface Methods
    • Ch13 Designing Robust Products and Processes: An Introduction
    • Ch14 Process Control, Forecasting, and Time Series: An Introduction
    • Ch15 Evolutionary Process Operation
  • Ch1
    • 考えること
    • データを取ること
    • データから考えること
    • そのループ
  • Ch2
    • (単)変量統計の基礎
  • Ch3
    • 2つのものを比較する
    • 推定する・検定する
  • Ch4
    • 複数のものを比較する
    • そのために必要になるデザイン上のこと
    • ブロックとランダマイズド・デザイン
    • ラテン方格と不完全ブロックデザイン
  • Ch5
    • 「要素のありなし」を複数要素についてデザインする
    • 2^k
  • Ch6
    • 「要素のありなし」を2^kにしないで、「はしょる」手立て
    • 2^{k-1}とか->拡張される
    • バランスよく「はしょる」ために、ダミー変数を使う
    • 「はしょっ」てあるだけに、「完璧」ではない
    • 要因のスクリーニングなどに用いればよし
    • 「はしょる」ときの原則は、「バランスよく」なので、それは、要因同士が「独立」と仮定しているということの裏返し
    • したがって、「独立」してそうなものについて、ダミー変数で組み合わせるのが良い(はず)
  • Ch7
    • 2^tに合わないときに、うまく「はしょっ」た組合せがほしい
    • 「バランスよい」ことは、原則
  • Ch8
    • Factorial designに鑑みて、データ変換のことを
  • Ch9
    • この本のデータの処理は、ばらつきが問題になる話だが、そもそもそのばらつきの起源について
  • Ch10
    • 回帰と最小二乗法にまつわること
  • Ch11
    • Ch12のResponse Surface Methods(RSM)を理解するための多次元幾何学的記載
  • Ch12
    • 「因子の働いている空間のどの部分」を調べたいのかを幾何学的にとらえよう
  • Ch13
    • 要因の関与の仕方を知ることで、頑健な工程を設計する
  • Ch14
    • 時系列の関与
    • 観察を時系列で行うと「モニタリング」
    • 時系列で行動を変化させれば「調節・調整」
  • Ch15
    • 日々の運用をしながら、情報を取得する(実験デザインではなく)

2011-07-29 AMMI Analysis

[][][][][]Additive Main Effects and Multiplicative Interaction Models

  • Association schemes -> 実験計画法 -> Rで実験計画法 -> agricolae パッケージ ときて、AMMI解析に出会う
  • AMMIの論文はこちら
    • 環境の違うところで実験をして、着目因子が、特性値にどのような影響力があるかをANOVAの枠組みで解析しようという手法

2011-07-27 実験計画法続き

[][][][]実験計画法メモ

  • 実験計画法入門byNAG
  • フィッシャーによる実験計画法
    • 実験計画法では研究を2つに分ける
      • 計画性を持ってデータを集まる段階
      • 集めたデータを解析する段階
  • 実験計画立案時のポイント
  • 集めたデータを解析する手法
    • 分散分析
    • 回帰分析
    • 多重比較

[][][][]実験計画法メモ2

  • 1元配置
# 1元配置
Ng<-5

Ns<-sample(5:10,Ng,replace=TRUE)
Fs<-c()
for(i in 1:Ng){
	Fs<-c(Fs,rep(i,Ns[i]))
}
Vs<-rnorm(sum(Ns))

D<-cbind(Fs,Vs)
print(D)
plot(Vs ~ Fs)
myfit <- lm(Vs ~ Fs)

anova(myfit)
  • 2元配置
    • 各水準の繰返し数によって「分散の分解」の成立関係が異なるので、対応が変わる
      • (1)各水準の繰り返し数が等しく,1 である場合
      • (2)各水準の繰返し数が等しく,2 以上である場合
      • (3)各水準の繰返し数が等しくないが,周辺度数に比例する場合
      • (4)各水準の繰返し数が等しくなく,周辺度数にも比例しない場合
# 2元配置
#(1) 各水準の繰り返し数が1
Ngs<-c(5,4)
Gs<-expand.grid(1:Ngs[1],1:Ngs[2])
Ng<-length(Gs[,1])
Ns<-sample(1,Ng,replace=TRUE)
Fs<-matrix(0,0,2)
for(i in 1:Ng){
	tmp<-matrix(c(rep(Gs[i,1],Ns[i]),rep(Gs[i,2],Ns[i])),ncol=2)
	Fs<-rbind(Fs,tmp)
}
Vs<-rnorm(sum(Ns))
D<-cbind(Vs,Fs)
print(D)
plot(Vs~Fs[,1]+Fs[,2])
myfit<-lm(Vs~Fs[,1]+Fs[,2])
anova(myfit)

#(2)(3)各水準の繰り返し数が周辺度数に比例する場合(すべて等しい場合を含む)

Ngs<-c(5,4)
Gs<-expand.grid(1:Ngs[1],1:Ngs[2])
Ng<-length(Gs[,1])
Ns<-sample(5,Ng,replace=TRUE)
Fs<-matrix(0,0,2)
for(i in 1:Ng){
	tmp<-matrix(c(rep(Gs[i,1],Ns[i]),rep(Gs[i,2],Ns[i])),ncol=2)
	Fs<-rbind(Fs,tmp)
}
Vs<-rnorm(sum(Ns))
D<-cbind(Vs,Fs)
print(D)
plot(Vs~Fs[,1]+Fs[,2]+Fs[,1]:Fs[,2])
myfit<-lm(Vs~Fs[,1]+Fs[,2]+Fs[,1]:Fs[,2])
anova(myfit)

#(4) 各水準の繰り返し数が周辺度数に比例しない場合
Ngs<-c(5,4)
Gs<-expand.grid(1:Ngs[1],1:Ngs[2])
Ng<-length(Gs[,1])
Ns<-sample(5:10,Ng,replace=TRUE)
Fs<-matrix(0,0,2)
for(i in 1:Ng){
	tmp<-matrix(c(rep(Gs[i,1],Ns[i]),rep(Gs[i,2],Ns[i])),ncol=2)
	Fs<-rbind(Fs,tmp)
}
Vs<-rnorm(sum(Ns))
D<-cbind(Vs,Fs)
print(D)
plot(Vs~Fs[,1]+Fs[,1]*Fs[,2])
plot(Vs~Fs[,2]+Fs[,1]*Fs[,2])
myfit1<-lm(Vs~Fs[,1]+Fs[,1]*Fs[,2])
myfit2<-lm(Vs~Fs[,2]+Fs[,1]*Fs[,2])
# myfit1 の Fs[,2] に関する結果を要因Fs[,2]の有意性に使う
anova(myfit1)
# myfit2 の Fs[,1] に関する結果を要因Fs[,1]の有意性に使う
anova(myfit2)
# 交互作用に関する出力はmyfit1もmyfit2も同じ

  • 多元配置
# 多元配置
#(1) 各水準の繰り返し数が1
k<-4
Ngs<-sample(3:7,k,replace=TRUE)
NgsList<-list()
for(i in 1:k){
	NgsList[[i]]<-1:Ngs[i]
}
#Gs<-expand.grid(1:Ngs[1],1:Ngs[2],1:Ngs[3],1:Ngs[4])
Gs<-expand.grid(NgsList)
Ng<-length(Gs[,1])
Ns<-sample(1,Ng,replace=TRUE)
Fs<-matrix(0,0,k)
for(i in 1:Ng){
	tmp1<-c()
	for(j in 1:k){
		tmp1<-c(tmp1,rep(Gs[i,j],Ns[i]))
	}
	tmp<-matrix(tmp1,ncol=k)
	Fs<-rbind(Fs,tmp)
}
Vs<-rnorm(sum(Ns))
D<-cbind(Vs,Fs)
print(D)
xnam <- paste("Fs[,", 1:k, "]",sep="")
(fmla <- as.formula(paste("Vs ~ ", paste(xnam, collapse= "+"))))


plot(fmla)
myfit<-lm(fmla)
anova(myfit)
  • 乱塊法
    • 2要因での乱塊法は各水準の繰り返し数がすべて1の場合の2元配置分散分析と同じになっている(ここここ)
  • ラテン方格法
    • 3元配置ANOVAで、「配置の仕方」が特別な(ラテン方格を使う)方法
    • ラテン方格の数え上げは難しく、単純な式や、漸化式もない(らしい)→こちら

2011-07-26 Association schemes2

ryamada222011-07-26

[][][]実験計画法

  • Association schemesは元々、実験計画法(Wiki)の「要因」のコントロールの配置の仕方から来ている話題
  • その教科書を2冊
入門 実験計画法

入門 実験計画法

[]Association Schemesの周辺

  • 日本語の記事が少ないけれど、こんな文章は参考になる

[][]Association schemes

  • こちらの続き
  • Latin squareが作るAssociation schemeの行列表現を作ってみる
  • ラテン方格は『n 行 n 列の表に n 個の異なる記号を、各記号が各行および各列に1回だけ現れるように並べたもの』
  • 今、n^2個の要素について、そのペアをとって関係を考える
    • 「自身」という関係
    • 「第一関係」という関係
      • ラテン方格において、行が一緒か、列が一緒か、入っている値が一緒だったら「第一関係」。ただし、「自身」という関係の場合を除く
    • 「第二関係」はそれ以外
  • 上記3つの関係として、それぞれの関係にある場合に1、ない場合に0となるような(n^2) ¥times (n^2)行列を3つ作る
# Latin square
# ラテン方格の文字の入り具合
X<-matrix(c("A","B","C","D","D","A","B","C","C","D","A","B","B","C","D","A"),4,4,byrow=TRUE)
# n^2要素に数値のIDをつける
Y<-matrix(1:16,4,4,byrow=TRUE)

M<-matrix(0,4^2,4^2)
# 自身という関係は対角成分がすべて1の行列
A0<-diag(rep(1,4^2))
# ラテン方格の行が一緒、列が一緒、ラテン方格の文字が一緒の関係を
# R, C, Lに入れる
R<-C<-L<-matrix(0,4^2,4^2)

for(i in 1:16){
	tmp1<-which(Y==i,arr.ind=TRUE)
	for(j in 1:16){
		tmp2<-which(Y==j,arr.ind=TRUE)
		if(tmp1[1]==tmp2[1]){
			R[i,j]<-1
		}
		if(tmp1[2]==tmp2[2]){
			C[i,j]<-1
		}
		if(X[tmp1]==X[tmp2]){
			L[i,j]<-1
		}
	}
}
# 3条件のいずれかに当てはまるものをA1で1にする
A1<-R+C+L-3*A0
# すべての要素が1の行列
J<-matrix(1,16,16)
# A2は「残り」
A2<-J-A1-A0

A1%*%J
A2%*%J
  • グラフ表示する

v<-16

t<-seq(from=0,to=1,length=v+1)*2*pi
t<-t[-length(t)]

plot(cos(t),sin(t))
As<-list()
As[[1]]<-A0
As[[2]]<-A1
As[[3]]<-A2

for(i in 1:length(As)){
	for(j in 1:length(As[[i]][,1])){
		for(k in 1:length(As[[i]][1,])){
			if(As[[i]][j,k]==1){
				segments(cos(t[j]),sin(t[j]),cos(t[k]),sin(t[k]),col=i,lwd=1)
			}
			
		}
	}
}