Site admin: motohiro PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 7.4.33
HTML convert time: 0.038 sec.
#contents
* 互いに相関しているサンプルの分散の差の検定 2006 01 06 [#w75a0586]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.187 - p.188
* 偏相関係数の計算方法 2006 01 06 [#w512ce2d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 191 - p.192
lm で，対象としたい変数を取り除いたモデルを update で生成すると
その変数の 偏 SSR (regression sum of squares) が求まる。
これを SST で割り，平方に開けば良い。
* SST (total sum of squares) の計算方法 2006 01 06 [#z45d71f9]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 192
lm(X ~ 1)
* Anova 表の見方，RSS 残差平方和 2006 01 06 [#xe0ae3d7]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 192
RSS は回帰によって説明されない変動，長谷川「多変量」 p.22
* [S] S-Plus で locale のエラー [#bb3c3fa8]
Warning: Locale not supported for XmbTextListToTextProperty
Warning: Cannot convert XmString to compound text
http://s.isac.co.jp/FAQ/faq-0032.html を参照
http://S.isac.co.jp/docs-4.5/machdep_3.htmlを参照
その他，myTurboAMD.html の S-plus インストールを参照
* [S] detach rm の順番 2006 01 09 [#qd4b3d6a]
まず detach し，次に rm する。逆に行なうと，
Problem: Object "factories" not found
Use traceback() to see the call stack
とエラー表示
* ネストの式 2006 01 11 [#b93e19cd]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.213
y ~ A/B
y ~ A * A:B
y ~ A + B %in % A
に等しい。
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.135 - 136
x
複数の要因からなるデータの nest は，ある要因の水準ごとに
別の要因の傾きを算出する。そして
A) すべての切片が 0 か
B) ある要因のある水準に対する別の要因の傾きがすべて 0 か
を検定する。
ネストモデルについては三中先生のサイト
http://cse.niaes.affrc.go.jp/minaka/R/NestedANOVA.html
を参照
* 切片を除く．その解釈 2006 01 11 [#t4e8647c]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.214
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p. 132 下, p.136
y ~ x - 1
は説明変数が連続量の場合，切片を 0 にし，全体平均と 全体 deviance を算出するが，
すべての説明変数がカテゴリデータの場合，各レベルの平均を算出する。
( -1 が無ければ，全平均，また平均からの差になる)。
* scale location plot 2006 01 14 [#n7eb0bef]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.237 -
回帰分析の結果を plot の引数に与えた場合の第三のプロットは，
分散が定数であるかを確かめるのに有効。
* カテゴリベクトルの順序づけ ordered 2006 01 14 [#f2988fef]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.251
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.51
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.130
* 変数削減の方法 anova と summary の違い 2006 01 18 [#ad4c556a]
Dalgaard, S-Plus, p.181, p.155
summary() で表示された P 値は，表示順に関係なく削除可能。
これに対して anova() の場合は，表示順序が重要。
また高次の相互作用をチェックする場合，ある単独ファクターが有意でなくとも，
それを含む高次作用が有意の場合は，その単独ファクターは残す。
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.269, p.270
* 適合値と分散の対応 2006 01 18 [#f519a964]
要因A と要因Bの各レベル語とに組合せにおける実験回数が少ない場合
多少，適合値に対して分散が減少する傾向があっても問題ない。
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 264
* [S] AIC の計算方法 2006 01 27 [#h08b644e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.316
mode.1<- lm(growth ~ 1)
AIC(model)
* NULL model はパラメータを grand mean 以外推測しない 2006 01 30 [#c522fe96]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.221, p.331
model<- aov(y ~ 1)
model<- lm(y ~ 1)
* Helmert contrasts vs. Treatment contrasts 2006 01 30 [#i751ed31]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.333, p.340
y = u + b1 x1 + b2 x2 + b3 x3 + b4 x4
というモデルが想定される場合，五つのパラメータを想定する必要はない。
treatment mean は
b1 を 0 と見なし，x1 の平均を u とし，他の b は，u との差と考える。
これに対して Helmert Contrasts は，
全体平均を最初のパラメータとし
最初と二つ目の平均の平均と，最初の平均との差を第二のパラメータ
最初と二つ目の平均の平均と，最初と二つ目さらに三つ目の平均の平均との差を第三のパラメータ
最初と二つ目さらに三つ目の平均の平均と，全体平均の差を第四のパラメータ
またコントラストを手作業で設定するのは [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.143
他にもここを参照
* 繰り返しのない二元配置分散分析 2006 01 31 [#af9d9e29]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 335
* Block, Plot, Split-plot の違い 2006 02 01 [#w63cb6bb]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.346
* tapply と list の組み合わせ 2006 02 01 [#fcee3a67]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.260, p.347
attach(splityield)
names(splityiled)
> names(splityield)
[1] "yield" "block" "irrigation" "density" "fertilizer"
だとして
tapply(yield, list(block, irrigation), sum)
等のように使う。
list については他に crawley, S-Plus, p.406, p.414
* Interaction が有意の時，main effect は削除できない 2006 02 02 [#p283a3dc]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.348
* Fixed Effects はそのファクターレベルに情報がある 2006 02 02 [#a49371ac]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.214 - 216, 350.
* [S] interaction.plot の使い方 2006 02 02 [#l37c66ab]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 350
* [S] model.table の使い方 2006 02 02 [#s2e21fb1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.253, p.271, p. 351
* Block と Plot, Split-Plot 2006 02 03 [#qc63bace]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 353
Block (= Plot) をさらに半分にして Split-Plot と呼ぶ。
* Treatment structure, Error structure 2006 02 03 [#efa123e7]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 354
* Predict と estimate の違い 2006 02 03 [#u93595c1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.361
また fitted と predict の違いは [[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.36
* Corrected, uncorrected Sum os Squares 2006 02 03 [#g3394f1c]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 227, p. 366
ある変数の値 に その変数の値 をかけた値の総和 uncorrected sums of squares
から
その変数の総和の二乗をその変数の数で割ったものを引くと corrected sums of squares
Corrected, uncorrected Sum os Squares の違いは，前者は変数の値からその平均を引いた偏差の二乗だが，後者は，変数の値そのものの自乗．
* 変数変換について 2006 02 07 [#a1daf86a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 394
* predict の使い方 2006 02 06 [#ra24fe2d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 405
* sample with replacement 2006 02 08 [#v2aabfa4]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.478
sample with replacement, 取り出した石をもとに戻す
sample without replacement もとに戻さない Hypergeometric
* ヒストグラムで整数値をバーの中心にする 2006 02 09 [#d601c0a0]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.489
break points を 0.5 間隔にする。
hist(mass, breaks = -0.5:16.5)
* Deviance (乖離度) について 2006 02 11 [#x276d817]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.29, p.122 に的確な説明
単独モデルの場合
total deviance と residual deviance は
それぞれ satured model との差を表している。
複数モデル比較の場合
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.30
「パラメータの少ない方のモデルがフィットしている」が帰無仮説
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.508, p. 516, p.526, p.539
regression 等ではモデルの当てはめ具合を SSE で測ったが，
GLM では，binomial な proportion データについて deviance を用いる。
また p.526
Deviance is equivalent to sums of squares in linear models
p.539 には
total deviance と residual deviance の違い。
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.58
に単独のモデルについての deviance と，
二つの入れ子になったモデル(どちらかがパラメータが多い)を比較する場合の deviance
について解説がある。
* GLM 小サンプルの weight を小さくする 2006 02 14 [#n8ef695a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.525
fail <- total - success
y <- cbind(success, fail)
model1 <- glm(y ~ explanatory,binomial)
* カテゴリ変数のレベル変更 (アルファベット順以外) 2006 02 14 [#hcbf091d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.251
ordered(target, levels = c("E","D","C","B","A"))
* binomial, exponential, poisson エラーについて(overdispersion 2006 02 14 [#jfcae797]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 526
Everit & Hothorn p.103 - 105.
dispersion parameter は １ にセットされる。
なぜなら binomial, exponential, poisson では分散は，データ直接計算されるのではなく，平均の関数として定義されているからである。
residual deviance が residual degrees of freedom よりも大きい場合は
overdispersion
この時，モデルの比較はデビアンスではなく，F検定量で行う ([[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.45)
dispersion パラメータの推定方法は faraway, p.47, p.60
* logistics analysis of covariance with proportion data 2006 02 15 [#p340ac2f]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.532
例えば三つの薬品のそれぞれについて，分量を変えながら，効果の割合を比較する。
* wilcox.testが有効ではないケース 2006 02 17 [#q7dce0b3]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.545
0,1 からのみなるようなデータの場合同順位が多くなる。
* G test，chisq とは別の方法による独立性の検定 2006 02 18 [#g9ebb92f]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.548, p.577
glm.crss<- glm(c(22,61,69) ~1, family=poisson)
1- pchisq( glm.crss[13]$deviance, glm.crss[7]$df.residual)
あるいは
glm.crss<- glm( c(15,10,10,15) ~ 1, poisson)
summary( glm.crss)
1 - pchisq(glm.crss$deviance, glm.crss$df.residual)
[1] 0.5695988]
これは
pchisq(glm.crss$deviance, glm.crss$df.residual, lower = FALSE)
と同じこと
* nuisance 変数 2006 02 18 [#gb4dd871]
crawley, p.551
* 高次の交互作用がある場合，それより低次の要因は省けない 2006 02 18 [#i6c86640]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.553
* モデル比較を多数くり返す場合の p 値 2006 02 18 [#n9a361c2]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.553
"With this many comparisons, I would always work at p = 0.01 for significance,
rather than p = 0.05 "
* 二つの要因に有意差があると出力されても，それは Intercept との差 2006 02 22 [#l387de4b]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.580
> summary(p.580.model)
Call: glm(formula = cbind(relief, patients - relief) ~ treat, family = binomial)
Coefficients:
Value Std. Error t value
(Intercept) -1.067841 0.2471428 -4.320744
treatB 1.959839 0.3427433 5.718094
treatC 2.468734 0.3666004 6.734128 これは Intercept と差
* GAM (Generalised Additive Models 2006 02 24 [#x27432da]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.602
y ~ s(x) + s(y) + s(z)
y ~ s + lo(w) + z
s は s smoother , lo は lo smoother を表す。
* Trend がある 2006 02 28 [#v681bb1e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.640
回帰分析をして，傾きに有意差があれば，そこにトレンドが認められる。
* プロット描画領域にテキストを加える 2006 03 02 [#wa2f9ccc]
text 関数を使う。例えば
text(grexp.obs$N[nrow(grexp.obs) ] , grexp.obs$K[nrow(grexp.obs) ], labels = "grexp")
これはデータフレーム grexp.obs の N 列の最大値を X 軸に， grexp.obs の K 列の最大値を Y 軸にとり，その交点にテキスト grexp を表示する。
~/research/statistics/baayen.R を参照。
* スケールの違うプロットを重ねる 2006 03 03 [#p0db652a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.656
par(mar = c(5,4,4,5) + 0.1) としてマージンを広げ
tsplot(cinnabar, ylab = "cinnabar") 最初のグラフを描き
par(new = T, xaxs = "d") として新しいプロットが先のプロットを上書きしないように設定
* 空のリストを作成する方法 2006 03 03 [#n3bb2f43]
Ligges p. 103
X <- vector(mode = "list", length=10)
* リストに対する演算 2006 03 03 [#w35e23a7]
X <- vector(mode = "list", length=5) # 5個の要素からなるリストを作成し
lapply(lapply(compare.K, "[", , 2 )
lapply(lapply(compare.K, "[", , "N" )
# リストの要素(データフレーム)のすべてから，すべての行の 2 列目(あるいはラベル名が N の列)を取り出す
lapply(lapply(compare.K, "[", , 2 ) , "max")
# リストから上記の条件で抽出した各要素について，その最大値を取り出す
max(as.numeric (lapply(lapply(compare.K, "[", , 1 ) , "max") ))
# 最終的に，リスト全要素の 2 列目の最大値を取り出す
* sapply の使い方 2006 03 04 [#k8e3c932]
Z.list の各要素はベクトルだとする。
unlist(Z.list) でベクトルが返される。
ちなみに sapply(Z.list, "[",1)
は最初の要素をベクトルとして取り出す。
sapply(Z.list, "[") はリストの各要素(ベクトル)を列とした行列を生成する。
Z.data<- data.frame(sapply(Z.list, "["))
とすればデータフレームに変更可能
* paste, parse, expression, eval によるデータフレームの初期化 2006 03 04 [#x4f04b03]
x<- "col"
y<- 1:5
z<- "th = 0"
xyz<- paste(x,y,z,sep="")
# > xyz
# [1] "col1th = 0" "col2th = 0" "col3th = 0" "col4th = 0" "col5th = 0"
xyz<- paste(xyz, collapse=",")
# > xyz
# [1] "col1th = 0,col2th = 0,col3th = 0,col4th = 0,col5th = 0"
d<- "data.frame("
d1<- ")"
d2<- paste(d,xyz,d1)
# > d2
# [1] "data.frame( col1th = 0,col2th = 0,col3th = 0,col4th = 0,col5th = 0 )"
parse(text =d2)
# > parse(text =d2)
# expression(data.frame(col1th = 0, col2th = 0, col3th = 0, col4th = 0,
# col5th = 0))
Z<- eval(parse(text =d2) )
# > Z<- eval(parse(text =d2) )
# > Z
# col1th col2th col3th col4th col5th
# 1 0 0 0 0 0
* 文字列の処理，文字列を式に変更 2006 03 04 [#s9534d8b]
Ligges, p.54, p.91
例えばデータフレームの列名をベクトルに変換するには
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.682 の例から
pig<- read.table("/home/ishida/source/statistics/S/crawley/data/pig.txt",header=T)
attach( pig )
names( pig )
deparse( names( pig ))
# "c(\"Pig\", \"t1\", \"t2\", \"t3\", \"t4\", \"t5\", \"t6\", \"t7\", \"t8\", \"t9\")"
parse(text = deparse( names( pig )))
#expression(c("Pig", "t1", "t2", "t3", "t4", "t5", "t6", "t7", "t8", "t9"))
eval(parse(text = deparse( names( pig ))))
#[1] "Pig" "t1" "t2" "t3" "t4" "t5" "t6" "t7" "t8" "t9"
parse(text = eval(parse(text = deparse( names( pig )))))
#expression(Pig, t1, t2, t3, t4, t5, t6, t7, t8, t9)
よって次のように行なう
pig.colnames < - eval(parse(text = deparse( names( pig ))))
* デ−タフレ−ムやリストをベクトルに変換 2006 03 06 [#n4d26bf3]
x1<- data.frame(y1 = c(1,2,3), y2 = c(11,22,33))
x2<- data.frame(y1 = c(4,5,6), y2 = c(44,55,66))
x3<- data.frame(y1 = c(7,8,9), y2 = c(77,88,99))
X123<- list(x1,x2,x3)
リスト X123 の各要素はベクトルだとする。
unlist(X123) でベクトルが返される。
# データフレームの各列をつなげる
unlist(x1)
あるいは
attach(x1)
z<- NULL
for(i in 1:2){
z<- c(z, get(paste("y", i ,sep="")))
}
# 列数が少ないなら
c(y1,y2) #でも OK
* rep()関数の使い方 2006 03 06 [#w4f9f6e7]
> rep(c("Carroll","Dickens", "Eliot", "Bronte", "Stevenson", "Gaskell"), rep(3,6))
あるいは
> rep(c("Carroll","Dickens", "Eliot", "Bronte", "Stevenson", "Gaskell"),each = 3)
[1] "Carroll" "Carroll" "Carroll" "Dickens" "Dickens" "Dickens"
[7] "Eliot" "Eliot" "Eliot" "Bronte" "Bronte" "Bronte"
[13] "Stevenson" "Stevenson" "Stevenson" "Gaskell" "Gaskell" "Gaskell"
* 文字列操作 2006 03 06 [#m580f36e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.37
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.305
Ligges p.54
例えば strtrim(labs,2) は文字列ベクトル labs のそれぞれの要素の最初の二文字を取り出す。
* アルファベットの略語 2006 03 07 [#w7fcf51b]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.38
abbreviate(str)
* 対数尤度比を算出 ，松原「統計学の考え方」 2006 03 10 [#z94fe472]
松原 「統計学の考え方」, p82 - p.84
> ftable(p.73)
death Y N
suspect victim
Wh wH 19 132
bL 0 9
Bl wH 11 52
bL 6 97
m2 iterations: deviation 8.881784e-16
$lrt
[1] 0
$pearson
[1] 0
$df
[1] 0
$margin
$margin[[1]]
[1] 1
$margin[[2]]
[1] 2
$margin[[3]]
[1] 1 2
> loglin(matu, c(1,2), param =T)
2 iterations: deviation 0
$lrt
[1] 0.1658080
$pearson
[1] 0.1690821
$df
[1] 1
あるいは
library(MASS)
matu<- array(data=c(4,6,20,40),
dim=c(2,2), dimnames = list(Row =c("a1","a2"), Col = c("b1", "b2")) )
loglm(~Row + Col, data = matu)
x<- array(data=c(15, 10, 10, 15), dim=c(2,2))
# = x<- matrix(c(4,6,20,40), nrow = 2)
kekka<- loglin(x, c(1,2), param =T)
1 -pchisq(kekka$lrt,kekka$df)# 対数尤度比 G による検定結果
1 -pchisq(kekka$pearson,kekka$df)# カイ二乗検定による検定結果
* R での時系列分析 2006 03 10 [#f9c56cde]
R-wiki
R の基本パッケージ base, stats 中の時系列オブジェクトの簡易解説
時系列データに関する最も基本的な統計量は
自己共分散(auto-covariance)と
自己相関係数(auto-correlation)である。
R の関連する関数は acf(), pacf(), ccf() である。
関数 acf() は時系列オブジェクトの自己共分散と自己相関係数を計算し、既定でそれをプロットする。
関数 pacf() は 偏自己相関係数(partial autocorrelations)を計算する。
関数 ccf() は二つの一次元時系列間のクロス相関係数(cross-correlation)と
クロス共分散(cross-cavarinace)を計算する。
acf(x, lag.max = NULL,
type = c("correlation", "covariance", "partial"),
plot = TRUE, na.action = na.fail, demean = TRUE, ...)
pacf(x, lag.max = NULL, plot = TRUE, na.action = na.fail, ...)
ccf(x, y, lag.max = NULL, type = c("correlation", "covariance"),
plot = TRUE, na.action = na.fail, ...)
* 回帰分析の summary の解釈 2006 03 14 [#c2d6ebc2]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.134
> summary(prestige.mod.3)
Call:
lm(formula = prestige ~ income + education + type + income:type +
education:type)
Residuals:
Min 1Q Median 3Q Max
-13.462 -4.225 1.346 3.826 19.631
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.276e+00 7.057e+00 0.323 0.7478
income 3.522e-03 5.563e-04 6.332 9.62e-09 ***
education 1.713e+00 9.572e-01 1.790 0.0769 .
typeprof 1.535e+01 1.372e+01 1.119 0.2660
typewc -3.354e+01 1.765e+01 -1.900 0.0607 .
income:typeprof -2.903e-03 5.989e-04 -4.847 5.28e-06 ***
income:typewc -2.072e-03 8.940e-04 -2.318 0.0228 *
education:typeprof 1.388e+00 1.289e+00 1.077 0.2844
education:typewc 4.291e+00 1.757e+00 2.442 0.0166 *
例えば，ここでは，Intercept は type要因の bc 水準で，これがベース
income, education は type要因の bc 水準の時の傾き
typeprf, typewc は切片との差
income:typeprof は，type要因の prof 水準の時の傾き
Dalgaard, p.178
Call:
lm(formula = log10(diameter) ~ log10(conc) * glucose, data = hellung)
Residuals:
Min 1Q Median 3Q Max
-2.672e-02 -4.888e-03 5.598e-05 3.767e-03 1.761e-02
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.627926 0.033754 48.230 < 2e-16 ***
log10(conc) -0.046716 0.006846 -6.823 1.51e-08 ***
glucose 0.003418 0.023695 0.144 0.886
log10(conc):glucose -0.006480 0.004821 -1.344 0.185
まず最初の二つは glucose があるデータの切片と傾き，次の二つはない場合の切片と傾き
二つの切片の差は 0.003418 , 傾きの差は -0.006480
差は二つとも有意でない(0.88, 0.185)
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.677
options(contrasts=c("contr.treatment","contr.poly"))
Call: lm(formula = Weight ~ Age * Genotype)
Residuals:
Min 1Q Median 3Q Max
-0.7787 -0.3986 -0.01139 0.4055 0.8445
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 7.5407 0.3853 19.5724 0.0000
Age 0.2931 0.1162 2.5230 0.0150
GenotypeCloneB 1.0401 0.5449 1.9090 0.0623
GenotypeCloneC -0.9487 0.5449 -1.7411 0.0881
GenotypeCloneD 0.5883 0.5449 1.0797 0.2857
GenotypeCloneE -0.9587 0.5449 -1.7596 0.0849
GenotypeCloneF 1.5693 0.5449 2.8802 0.0059
AgeGenotypeCloneB -0.0241 0.1643 -0.1468 0.8839
AgeGenotypeCloneC -0.0316 0.1643 -0.1926 0.8481
AgeGenotypeCloneD 0.0786 0.1643 0.4782 0.6347
AgeGenotypeCloneE 0.0278 0.1643 0.1690 0.8665
AgeGenotypeCloneF -0.0116 0.1643 -0.0705 0.9441
ここで (Intercept) はGenotypeCloneA を基準とした ベースとなる切片
Age はベースとなる傾き
AgeGenotypeClone* は交互作用で，ベースとなる傾きに対する差
ここも参照
また多重回帰分析の場合は，切片以外の項はその変数に対する傾き
* 文字列ベクトルを数値ベクトルにかえる, 要因をまとめる 2006 03 14 [#v224940f]
as.numeric(as.factor(string.vector))
要因をまとめる crawley, S-Plus, p.301
newgen<- factor(1+(Genotype=="CloneB")+(Genotype=="CloneD")
+ 2*(Genotype=="CloneC")+2*(Genotype=="CloneE")
+ 3*(Genotype=="CloneF"))
lme における anova の解釈 2006 03 17
crawley, S-Plus, p.679, p.690 - 691
anova(k.model1, k.model2 )
# crawley, p.679 複雑なモデル k.model2 が AIC が小さくて better
# 説明力( ## Model df AIC BIC logLik Test L.Ratio p-value
#k.model1 1 3 2525.512 2537.162 -1259.756
#k.model2 2 4 1932.097 1947.630 -962.048 1 vs 2 595.4153
しかし AIC が低くても，p-value が大きければ，説明力はない
* [S] S-Plus での画像出力 2006 03 17 [#a03d584c]
http://www.math.aau.dk/~dethlef/Links/latex_figures.html
There are two possibilities for directly producing a PostScript version
of an Splus graphic:
use the menu in the graphics window, with the print command set to "echo",
so the graph is saved to a file (ps.out.001.ps etc.)
instead of being sent to the printer.
Alternatively, from a program use something like:
postscript(file="foo.ps", horiz=F, onefile=F, print.it=F)
plot...
dev.off()
unix("ps2epsi /home/ishida/mysplus/foo.ps /home/ishida/mysplus/foo.eps")
The output will be written to the specified file.
It may be useful to process the PostScript file
with the ps2epsi program (supplied with ghostscript)
before including it into latex, since Splus tends to leave large margins:
ps2epsi foo.ps foo.eps
This also adds a bitmap version of the image,
which is apparently used by certain Macintosh software to display it
on the screen.
* [S] unix への操作 2006 03 17 [#e664fa12]
http://www.uni-muenster.de/ZIV/Mitarbeiter/BennoSueselbeck/s-html/helpfiles/unix.html
# create a character vector from the contents of a file,
# one element of the vector per line of the file
unix("cat myfile")
unix("cat myfile",output=F) # same as
!cat myfile
unix("more", c("line 1", "here is line 2", "and line 3"), F)
# C shell commands will work if S_SHELL is set to "/bin/csh"
# (You can set S_SHELL explicitly, or let it default to SHELL or /bin/sh)
!alias # works if S_SHELL is /bin/csh
# You can use unix.shell in place of unix to run csh commands
unix.shell(command = "alias", shell = "/bin/csh", out=F)
# You can make a convenient function to run csh commands
"csh"<-
function(cmd, ...)
{
unix.shell(command = cmd, shell = "/bin/csh", ...)
}
* [S] 余白の大きさを変える 2006 03 17 [#gaa6256e]
http://www.msi.co.jp/splus/support/salon/mcourse/mc12.html
par(mar=c(2, 3, 2, 2))
# 余白の大きさを，行(フォント)の高さを単位に変更
などです。
これらパラメータを設定されると、次に par で同じパラメータが呼び出されるまで、
あるいは、グラフ・ウィンドウが閉じられるまで継続することを忘れないでください。
また S-plus 日本語訳の p.81
ただし lattice パッケージの trellis.device については
ligges, p.167
* [S] Java アプレットでのオンラインヘルプ 2006 03 18 [#gb7baad8]
S-plus のコマンド状で
help.start()
Starting up Java help system. This may take a minute.
終えるには
help.off()
* Regression constant (切片)を除く指定と，その意味 2006 03 18 [#y66b9264]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.453 - p.454, R: p.109
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.132, p.136
切片(総平均)を除くにはモデル構造式に -1 を含めるが，その意味は
ある要因のそれぞれのグループ(水準)に，異なった切片をフィットさせることである。
* 要因を数値に変える unclass 2006 03 18 [#vcc6997d]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.127
* 様々な Anova I, II, III 2006 03 20 [#ra068f7c]
http://www.aichi-gakuin.ac.jp/~chino/anova/chapter1/sec1-3-7.html
他に [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.135
SAS では、glm プロシジャで、分析に際して４つの 平方和が用意されており、
Type I 平方和、 Type II 平方和、Type III 平方和、Type IV 平方和と呼ばれる。
要因が ２つ以上あるデザインの場合、これらは場合によって異なるので、
その選択には 注意を要する。
まず、Type I 平方和は、逐次平方和 (sequential sums of squares) とも呼ばれ
(Draper & Smith, 1981)、複数の要因を並べた順に 追加していく時のモデル平方和
の増加を評価する。したがって、Type I 平方和は一 般的には、要因の投入順序に
より異なった値を取る。ただし、要因相互が直交して いる場合
（主として釣り合い型デザインの場合）には、この平方和は投入順序に依 存しない。
Type II 平方和は、当該要因と交互作用を持つような要因や同じく当該要因と
交絡するような要因以外のすべての要因の影響を差し引いた平方和で、
偏平方和 (partial sums of squares) とも呼ばれ る (Draper & Smith, 1981)。
したがって、 Type II 平方和は、要因の投入順序に は依存しない (SAS, 1990)。
Searle (1987) によれば、Type III 平方和は、Σ 制約付 モデルの
平方和 (SS for Σ-restricted models) である。一方、Type IV 平 方和は、
仮説平方和 (hypothesis SS) であり、SAS glm ルーチン自身により決定される
仮説検定のための平方和である。
Type III 平方和も Type IV 平方和も、偏平方和と呼ばれることがある (SAS, 1990)。
デザインの中に全く欠測値がない場合、両者は一致する。
釣り合い型デザインでは、４つの平方和は、すべて等しい。また、交互作 用
のない場合には、Type II 平方和、Type III 平方和、Type IV 平方和は一致する
(Searle, 1987)。また、非釣り合い型デザインでも、交互作用項を含まない主効果
のみのデザインであれば、Type II 平方和とType III 平方和は一致するが、
交互作用項を含んだデザインでは、両者は異なる（竹内監修、高橋ら、1990)。
このような場合、高橋らはType II 平方和を用いることを薦めている。
* Jittering で重複点をずらしてプロットする 2006 03 20 [#ub4b7f94]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.139
points(jitter(Fcat[partner.status == 'low']), conformity[partner.status == 'low'], pch = 'L')
points(jitter(Fcat[partner.status == 'high']), conformity[partner.status == 'high'], pch = 'H')
* 要因Aと他要因 BCD全体との差，BCD 内の差を一度に検定 2006 03 20 [#c4f205d7]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.142 - 143.
* 欠損値の扱い 2006 03 21 [#s6421c0b]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.151
na.action(引数) で設定するが，デフォルトは
R は na.omit (欠損値を含まないデータで解析)
S は na.fail (エラーを出して解析を中止)
S では，計算の際に引数として na.action = na.omit, na.action = na.exclude
が必要。この二つの違いについては [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.151
* offset の意味 2006 03 21 [#h6f113a3]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.153
R のみでのオプション
線形モデルでは，目的変数からその変数を引くのに等しい。
faraway, GLM, p.63
rate mode で使う
* 対話的にテキストを挿入 locator 2006 03 22 [#if668996]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.165
text(locator(2), c("Close", "One-Sided"))
* contr.poly の出力の L, Q について 2006 03 22 [#ua736596]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.130
* car ライブラリの expand.grid 2006 03 23 [#n73552d4]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.170, p.181
データフレームの行数に合わせて，各列の変数を繰り返す
* [S] テーブルをデータフレームへ変換 2006 03 24 [#r2312161]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.184
cbind(expand.grid(dimnames(table)))
Freq = as.vector(table)
* logit と probit の違い 2006 03 28 [#sb64a734]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p. 27,
http://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%83%83%E3%83%88
Logit(p) = log(p/(1-p)) = log(p) - log(1 - p)
http://www.geocities.com/hikachu/2BAStatar.html
Logit&Probit
従属変数が0-1の値をとる場合の回帰分析です．
そのまま通常の回帰分析をすると予測値が0-1の範囲を飛び出てしまったりしてまずいので，
それを避けるために0-1の従属変数に特殊な変換を施して[-∞,+∞]になるようにしてやります．
こうしておいて好きなだけ推定してから逆の変換を施します．
そうすれば予測値は見事[0,1]に収まるというわけです．
この変換にロジスティック関数を用いるのがLogit，累積正規分布関数を用いるのがProbitです．
本当はNormitの方が相応しい名前かもしれませんが…．
LogitとProbit，どういう基準で使い分けるべきか?というのが初心者共通の悩みですが，
これはどちらでもいいです．どちらを用いても違いはありません．同じです．
若干の違いがあるとすれば 1.レアなイベント 2.拡張版の相性 ぐらいです．
1について．累積正規分布関数よりロジスティック関数の方がtailがheavyなので
自分の関心がレアなイベントにあるときはProbitよりLogitの方がよいとか．
もっとtailがheavyな関数を使ったものとしてCauchyというのもあります．
2について付け加えると，例えば，multiple probit, heteroscedastic logitは
技術的に困難だそうです．
* drop1 は説明変数単独の有意度を測る 2006 03 28 [#ue93c526]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.33
full モデルとの差を測る
* 信頼区間を求める confint 2006 03 28 [#z9e323cc]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.34
* R with 関数 2006 03 30 [#u4cb5e2d]
以下はRjpwikiより
データを引数 data = hoge と指定できない場合などに便利
例えば interaction.plot を attach されていないデータ troutegg に適用する場合
with(troutegg, interaction.plot(period, location, elogits))
を実行する。
interaction.plot(period, location, elogits, data = troutegg)
* ケースとコントロール 2006 03 31 [#p2df496e]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.48
* バートレットの多群の等分散性 bartlett.test 2006 04 01 [#zd61ecee]
永田靖著『統計的多重比較法の基礎』サイエンティスト社
bartlett.test(p.37$x ~ p.37$name)
* 群ごとに分散などを求める 2006 04 01 [#y7594d12]
nagata.aov.R 永田靖著『統計的多重比較法の基礎』サイエンティスト社
tapply(p.37$x, list(p.37$name), var)
* Dunnett の方法による多重比較 2006 04 01 [#sd223ce8]
これは対象群を中心に，それと処理群の平均値の差を調べる
永田 靖 p.40
p.43<- data.frame(x = c(7,9,8,6,9,8,11,10,8,8,
8,9,10,8,9,9,10,12,
11,12,12,10,11,13,9,10,
13,12,12,11,14,12,11,10),
name = (rep(c("A1","A2","A3","A4"), c(10,8,8,8)))
)
mulcomp ライブラリを使う
library(multcomp)
# p.43 のデータ形式のままだと A1 を引く形になる。
# そこでコントラストを変える
c.m<- matrix(c(1,-1,0,0,
1,0,-1,0,
1,0,0,-1), nrow = 3, byrow = T)
tapply(p.43$x, list(p.43$name), mean)
# 8.400 - 9.375 # -0.975
summary(simint(x ~ name, data = p.43, cmatrix = c.m, alternative = "less"))
* ロジット とは 2006 04 03 [#uc13557e]
http://www.clg.niigata-u.ac.jp/~takagi/basic90.html
対数オッズは，通常「ロジット logit」と呼ばれる．
log[ｐ／(1‐ｐ)]＝βｘ＋α
式を書き直すと，
ｐ＝[1＋exp{-(βｘ＋α)}]-1
と表わすことができる．
２群のロジットの差をとると，
logit(ｐ1)‐logit(ｐ2)＝log[ｐ1(1-ｐ2)/[ｐ2(1-ｐ1)]＝β
となる．
これは対数オッズ比log(ψ)がβであることを示している．
すなわち，ψ＝exp(β) である．
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.27, p.35
* /home/research/statistics 直下の eps ファイルの移動 2006 04 04 [#s51781bd]
/home/research/statics/EPS
を作成
* アルファベット順以外のカテゴリ分け関数 gl 2006 04 05 [#hcccef5c]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.70, p.101
カテゴリの基準をアルファベット順以外で構成する時などに使う
y<- c(320, 14, 80, 36)
particle<- gl(2, 1, 2*2, labels = c("no","yes"))
quality<- gl(2, 2, labels = c("good", "bad"))
wafer<- data.frame(y, particle, quality)
* GLM.summary の見方 2006 04 05 [#ae2f5560]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p. 70
> wafer
y particle quality
1 320 no good
2 14 yes good
3 80 no bad
4 36 yes bad
Call:
glm(formula = y ~ particle + quality, family = poisson)
Deviance Residuals:
1 2 3 4
1.324 -4.350 -2.370 5.266
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.6934 0.0572 99.535 <2e-16 ***
particleyes -2.0794 0.1500 -13.863 <2e-16 ***
qualitybad -1.0575 0.1078 -9.813 <2e-16 ***
> predict(modl)
1 2 3 4
5.693358 3.613916 4.635807 2.556366
> model.matrix(modl)
(Intercept) particleyes qualitybad
1 1 0 0
2 1 1 0
3 1 0 1
4 1 1 1
attr(,"assign")
[1] 0 1 2
attr(,"contrasts")
attr(,"contrasts")$particle
[1] "contr.treatment"
attr(,"contrasts")$quality
[1] "contr.treatment"
解釈は以下のように，ここで bad ではなく good が基準となっているのは gl 関数による
particle no, quality good 5.693358
particle yes, quality good 5.6934 - 2.0794(x1)
// この行に対応する Pr(z)は<2e-16 は particle no との差が有意なのを表す
particle no, quality bad 5.6934 - 1.0575(x2) /
/ この行に対応する Pr(z)は<2e-16 は quality good との差が有意なのを表す
particle yes, quality bad 5.6934 - 2.0794(x1) - 1.0575(x2)
> predict(modl,type = "response")
1 2 3 4
296.88889 37.11111 103.11111 12.88889
---
* outer 同時確率を出す 2006 04 05 [#a3dc4345]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.72
> wafer
y particle quality
1 320 no good
2 14 yes good
3 80 no bad
4 36 yes bad
> outer(qp,pp)
particle
quality no yes
good 0.6597531 0.08246914
bad 0.2291358 0.02864198
* deviance を Pearson X^2 で測る, すなわち chisq の別の方法 2006 04 05 [#v6fe92fe]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.73, p.76
> wafer
y particle quality
1 320 no good
2 14 yes good
3 80 no bad
4 36 yes bad
xtabs を使う．
(ov<- xtabs(y ~ quality + particle))
prop.test(ov)
* 特異値分解 2006 04 05 [#ha9a4eb1]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.77
関数 svd の対象は residuals
* plot 関数への asp 引数 2006 04 05 [#cef1e350]
http://www.ma.hw.ac.uk/ams/Rhelp/library/graphics/html/plot.window.html
Note that if asp is a finite positive value then the window is set up
so that one data unit in the x direction is equal in length to asp * one
data unit in the y direction.
* pchisq への lower = F 引数 2006 04 06 [#f890cd5f]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.80
pchisq(deviance(mods), df.residual(mods), lower =F)
# [1] 0.003762852
pchisq(deviance(mods), df.residual(mods))
# [1] 0.9962371
* margin.table の使い方 2006 04 06 [#f8f840f6]
< ct
left
right best second third worst
best 1520 266 124 66
second 234 1512 432 78
third 117 362 1772 205
worst 36 82 179 492
margin.table(ct, 1)
right
best second third worst
1976 2256 2456 789
* テーブルをデータフレームに変換 2006 04 06 [#ed1e2358]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.89
data(nes96)
# xtabs 使う
xtabs(~ PID + educ, nes96)
(partyed as.data.frame.table(xtabs( ~ PID + educ, nes96)))
* xtabs の使い方 2006 04 06 [#y76b4294]
xtabs(target.values ~ categorical.var1 + categorical.var2, data=data.name)
あるいは
> xtabs( ~ Q1 + Q2, data = mydata[,2:3])
Q2
Q1 1 2 3 4 5
1 3 2 3 2 0
2 0 0 2 0 0
3 0 0 1 1 1
4 0 0 0 0 1
* cut 関数の使い方 2006 04 06 [#fc4e5ea9]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.98
cutinc さらに ここ も参照
* log による summary 出力を antilog (exp) 化して解釈する 2006 04 07 [#u99d8e17]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.105, p. 34 に簡易的な見方
> summary(binmodw)
Call:
glm(formula = cbind(CNS, NoCNS) ~ Water + Work, family = binomial,
data = cns)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.4325803 0.0897889 -49.367 < 2e-16 ***
Water -0.0032644 0.0009684 -3.371 0.000749 ***
WorkNonManual -0.3390577 0.0970943 -3.492 0.000479 ***
ここで WorkNonManual の Coefficient -0.3390577 は次のように解釈できる。
log ベースで基準より -0.339 少ないとは，antilog に戻すと，
exp(-0.3390577) 0.7124413
であり，これは
基準 * 0.7124413, すなわち基準より約 29%少ないと解釈できる。
* P-value と null distribution 2006 04 08 [#caf4b323]
P-value と null distribution
http://hawaii.aist-nara.ac.jp/~shige-o/cgi-bin/wiki/wiki.cgi?QValue#i3
帰無仮説どおりのランダムサンプリングによって得られたデータの分布を、
null distribution と呼ぶ。
例えば、 T-検定を使用するような状況で t-統計量の null distribution は t-分布に従うし、
chi^2-検定を使用するような状況で 統計量の null distribution は chi^2-分布に従う。
ところが、どんな検定でも検定統計量のP-valueを計算すると、 そのnull distribution は
[1,0]の一様分布となる。 逆に、P-value はそうなるように定義されていると解釈してもよい。
つまり、 「P-value とは、その値が小さければ小さいほど仮説の有意性が高い
と解釈できるような指標であって、帰無仮説(null hypothesis)が真であるときに、
その値の出現頻度が[1,0]で一様分布するもの」である。さらに他の言葉で言えば、
「P-value は有意性指標を正規化したもの」である。
どんな帰無仮説に基いてどんな検定統計量を用いるどんな検定でも、
P-value が同じである仮説同士は同じぐらい有意 であると言うことができる。
* Wald test と，その問題点 2006 04 08 [#m78ef4d2]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.30, p.122 - 123
* データを中心化 (center) する 2006 04 15 [#f835bdd4]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.186
* iid 独立同型分布 2006 04 15 [#z3694c79]
iidは独立同型分布(independently and identically distributed)
* expand.grid 変数の値を変化させる(ただし他変数は固定) 2006 04 18 [#p85a707a]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.273
expand.grid(temp = seq(-3, 3, 0.1), ibh = 0, ibt = 0)
ibh, ibt の値は平均(標準化された平均値)に固定し，temp を[-3,3] の間で変化
* diagnostic.plot の見方 2006 04 19 [#c00ea4e5]
Wood, GAM, p.5 下 - p.6
* 任意のヒストグラムの描き方 2006 04 27 [#s35270d9]
例えば平均値が 64.68, 標準偏差が 13.97 のヒストグラムであれば
x<- c(64.58-(13.97*3) : 64.58+(13.97*3))
y<- dnorm(x, 64.68, 13.97)
plot(x,y, type="l")
lines(c(64.68, 64.68), c(max(y), min(y)))
lines(c(64.68 - 13.97, 64.68 - 13.97), c(0, .0169))
lines(c(64.68 + 13.97, 64.68 + 13.97), c(0, .0169))
text(66, -0.0005, "64.68")
lines(c(64.68 - 13.97, 64.68), c(.0169, .0169), lty = 3)
text(58, .0159, "13.97")
lines(c(min(x), max(x)), c(0,0))
* lm 関数出力の F 値 2006 04 28 [#qa233632]
Wood, p.30
推測されたパラメータが，切片しか含まないモデルから生成された確率。
* lm 関数出力の R^2 がマイナスのケース 2006 04 28 [#b3ea87f7]
Wood, p.33
モデルに不必要な説明変数が多く存在する
* A % in % B による条件判定とデータ抽出 2006 04 28 [#ucc36199]
Wood, p.35
sperm.comp1$m.vol<- sperm.comp2$m.vol[sperm.comp2$pair %in%
sperm.comp1$subject]
他に [[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.165
offset の意味 2006 04 28
Wood, p.45
モデルの含める際，そのパラメータが 1 と仮定する。
* layout 関数によるプロット座標の分割 2006 05 04 [#hb41b2cb]
attach(iris)
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE), c(2,1), c(1,2), TRUE)
layout.show(4)
for(i in 1:4){
boxplot(iris[, i] ~ Species, main = colnames(iris)[i])
}
Murrell, p.78
# layout を使うなら
layout(rbind(c(1,2),c(3,4)))
# layout.show(4)
# layout(matrix(c(2,0,1,3), 2, 2, byrow = TRUE), c(2,1), c(1,2), TRUE)
for(i in 1:4){
boxplot(iris[, i] ~ Species, notch = TRUE, main = colnames(iris)[i])
}
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.31 上
Wood, p.36
layout(matrix(c(2,0,1,3), 2, 2, byrow = TRUE), c(2,1), c(1,2), TRUE)
はプロットの横軸に関しては，左が全体の 2/3 右が 1/3 のサイズとし
縦軸に関しては，上が全体の 1/3 下が 2/3 とする設定である。
> matrix(c(2,0,1,3), 2, 2, byrow = TRUE)
[,1] [,2]
[1,] 2 0
[2,] 1 3
pairs と panel の使い方 2006 05 06
Hothorn, p.68
pairs(means[,-1], panel = function(x, y)
{text(x,y, abbreviate(levels(skulls$epoch)))})
* substitute() の使い方 2006 05 08 [#k083e5be]
substitute のヘルプより
f1<- function(x, y = x) { x<- x + 1; y }
s1<- function(x, y = substitute(x)) { x<- x + 1; y }
s2<- function(x, y) { if(missing(y)) y<- substitute(x); x<- x + 1; y }
a<- 10
f1(a)# 11
s1(a)# 11
s2(a)# a
> x<- 1
> substitute(x+1)
x + 1
> substitute(x+1, list(x = 2))
# x+1 の変数を指定された環境(ここではリスト)で置き換える
2 + 1
他に [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p. 831
* bubble plot の効果的な使い方 2006 05 12 [#uc43378a]
Hothorn, p.98
* multinom 関数の使い方 2006 05 17 [#udc862d4]
rmultinom(n, size, prob)
rmultinom(10, size = 15, prob=c(0.9,0.9))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 7 6 6 8 8 7 6 8 4 8
[2,] 8 9 9 7 7 8 9 7 11 7
10 は試行の回数,
15 は玉の最大個数
prob は箱の個数と，それぞれに入る玉の個数，ここでは箱は2個。
* glm 出力から Null deviance Residual deviance を抽出 2006 07 19 [#y628cb06]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.41
modl<- glm(cbind(dead, alive) ~ conc, family = binomial, bliss)
summary(modl)
modl$dev
modl$null
とする．
* ポアソン分布 2006 07 21 [#a53405f7]
工学のためのデータサイエンス, p.47 - 48
mase.47<- c(rep(0,55), rep(1,72), rep(2,40), rep(3,9), rep(4,5), rep(5,1))
mean(mase.47) #= (72 + 40 *2 + 3 * 9 + 4 *5 + 5) / 181
var(mase.47)
工業統計学, p.44
ppois(0,2)
ppois(1,2) - ppois(0,2)
* table 関数で作成したカウントのカテゴリ頻度 0 の調整 2006 07 25 [#y3c5dae0]
# 高頻度では、データ区間が飛び飛びになっているので、これを補正する
y<- min(buntyo):max(buntyo)
y<- 0:max(buntyo)
bun.df<- data.frame(cate = y, freq = c(rep(0, length(y))))
z<- 0
bun.df[1,] = c(0,0)
for(i in 1:nrow(bun.df)){
if(bun.orig.df[i - z,]$buntyo == bun.df[i,]$cate){
bun.df[i,]$freq<- bun.orig.df[bun.orig.df$buntyo == bun.df[i,]$cate, ]$Freq}
else{
z<- z + 1;
next;
}
}
##
* カイ二乗検定での期待度数 2006 07 27 [#oa50e416]
http://aoki2.si.gunma-u.ac.jp/lecture/Cross/warning.html
以下のいずれかに該当する分割表に対して χ2 検定を行ったり，
カイ二乗値から導き出される関連性の指標を用いる場合には注意が必要である
（Cochran, W. G.: Some methods for strengthening the common χ2 tests. Biometrics, 10, 417-451, 1954.）。
期待値が 1 未満の桝目が 1 つでもある。
期待値が 5 未満の桝目が全体の桝目の数の 20 % 以上ある。
このような場合には，
カテゴリーを併合する。
カテゴリーが併合できないような場合には，他の検定手法の適用の可能性を検討すべき
である（2 × 2 分割表の場合には「フィッシャーの正確確率検定」も参照）。
いずれも不可能な場合には，結果の解釈・適用において十分な検討が必要である。
教科書によっては 2 種類の χ2 検定が示されている。すなわち，「いくつかの群で，
ある変数の分布に差があるかどうかの検定」と，「2 変数が独立であるかどうかの検定」
である。両者は形式的には全く同じであるが，データの採取法により両者は区別される。
例えば，「男女の 2 群で食物嗜好の回答の分布が異なるか」を知りたいときに，
まえもって m 人の男子と n 人の女子を対象として調査することを決めてからデータを
集めた場合には男女で分布が異なるかどうかの検定になる。これに対して，N 人の対象者
について調査し，得られたデータを性別と回答別に集計した場合には，性別と食物嗜好が
独立であるかどうかの検定になる。疫学調査でいえば，前者は「後向き調査」や
「前向き調査」から得られる分割表の分析であり，後者は「断面調査」により得られる
分割表の分析である。
「二群の比率の差の検定（正規分布による検定）」は，2 × 2 分割表に対する
独立性の検定（χ2 検定）と等価な検定である。
検定結果が有意になるように（あるいは，逆に有意にならないように）
カテゴリーの分割方法（分割点）を変えたり，カテゴリーの併合を行ったりしては
いけない。
帰無仮説が棄却され「2 変数は独立ではない」となった場合にも，どのような関連があるかは
期待度数と観測度数を比較して考察しなければならない。
必ずしも線形の関連があるわけではない。
* ks.test を使った適合度の検定 2006 07 27 [#o50117f1]
x<- rnbinom(30, size = 5, prob = 0.5)
ks.test(x,"pnbinom", size = 5, prob = 0.5)
# 間瀬 p.51 下
n<- 4+9+16+13+9+7+5+4+3
x<- rnbinom(n, size = 10.59, prob = 0.7609)
ks.test(x, "pnbinom", size = 10.59, prob = 0.7609)
table(x)
ks.test(table(x), "pnbinom", size = 10.59, prob = 0.7609)
ks.test(observed, "pnbinom", size = my.k, prob = (my.k / (my.k + mean(x))))
ここで size は ２項分布のパラメータ k (= size)
prob は prob<- size/(size + mu) で計算される．
ks.test(rnorm(100), "pnorm")
ks.test(rpois(100, 5), "ppois", lambda = 5)
* glm を使った適合度の検定 2006 07 29 [#i4712bea]
null model を当てはめてみる
# データは負の二項分布に従っているか
x<- rnbinom(100, size = 2, prob=0.5)
x.glm.nb<- glm.nb(x ~ 1)
x.glm<- glm(x ~ 1, family = negative.binomial(1.57))
1 - pchisq( x.glm$deviance, x.glm$df.residual)
[1] 0.2081144
ポアソン分布ではどうか
x.glm.poi<- glm(x ~ 1, family = poisson)
1 - pchisq( x.glm.poi$deviance, x.glm.poi$df.residual)
[1] 4.327188e-09 # poisson による当てはめは悪い
他にも実験
# ポアソン分布を負の二項分布で当てはめる
x<- rpois(100, 3)
x.glm.nb<- glm.nb(x ~ 1)
summary( x.glm.nb)
1 - pchisq( x.glm.nb$deviance, x.glm.nb$df.residual)
# [1] 0.9046653
# 正規分布を負の二項分布で当てはめる
x<- rnorm(100, mean = 10, sd = 2)
x.glm.nb<- glm.nb(x ~ 1)
summary( x.glm.nb)
1 - pchisq( x.glm.nb$deviance, x.glm.nb$df.residual)
# [1] 1 当てはまってしまう
ここも
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.214, p.511, p.517, p.539 も
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.121 - 122.
このような記事もあった
また次のような記事もあった
From: Dan Kehler (kehler@mathstat.dal.ca)
Date: Sat 29 May 2004 - 02:54:48 EST
* Next message: Richard Valliant: "[R] dotchart questions"
* Previous message: Ted Harding: "RE: [R] Generate a sequence of random integer values"
Message-id:<Pine.GSO.3.96.1040528133631.18220B-100000@chase>
Using R 1.8.1, and the negative binomial glm implemented in MASS,
the default when using anova and a chi-square test is to divide the
deviance by the estimated dispersion. Using my UNIX version of S-plus (v
3.4), and the same MASS functions, the deviances are *not* divided by the
estimated dispersion.
Firstly, I'm wondering if anyone can enlighten about the correct procedure
(I thought the F-test was more appropriate when dispersion is estimated)?
Secondly, after a bit of muddling with the negative binomial pdf, I
concluded that, like for the Poisson, phi is actually 1. This result is
borne out by simulations. Is this correct?
# an example in R 1.81 with library(MASS)
y<-rnegbin(n=100,mu=1,theta=1)
x<-1:length(y) # 石田 追加，これは便宜上の説明変数を作ったに過ぎない
model<-glm(y~x,family=neg.bin(1))
summary(model)$dispersion
[1] 1.288926
anova(model,test='Chisq")
#...
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 99 102.038
x 1 0.185 98 101.853 0.705
# But the "real" chi-square probability is
1-pchisq(0.185,1)
[1] 0.6671111
1-pchisq(0.185/1.288926, 1) # 石田 追加
[1] 0.7047963 # 石田 追加
Thanks in advance,
Dan
これは，フルモデルにたいして，ヌルモデルを比べ，観測値の変動をヌルモデルだけで
どれだけ説明できるかを示したものである．
フルモデルは，本来，deviance がゼロのはずであるので，ヌルモデルとフルモデルの
デヴィアンスの差とは，ヌルモデルのデヴィアンスのことではないのか．
また
x<-1:length(y)
は単変量であって，これを使って回帰を行ってもフルモデルを設定したことには
ならないのではないか．そもそもフルモデルは自由度 ０ のはずである．
その意味で，ヌルモデルそのものがある分布に適しているかどうかを調べるには
先の方法でも十分かも知れない．
下の場合，もともと負の二項分布にのみ従った乱数で，他に変動要因が無い．
尤も，この場合，誤差は負の二項分布の誤差に従っているということになるのか?
y<- 1:length(x)
x.glm.nb2<- glm.nb(x ~ y)
anova(x.glm.nb2, test = "Chi")
これは，以下に等しい
pchisq( x.glm.nb$null.deviance -x.glm.nb$deviance ,
x.glm.nb$df.null - x.glm.nb$df.residual, lower = F)
anova(x.glm.nb, x.glm.nb2, test = "Chi")
以下の数値とは異なる数値が出る
x.glm.nb2$null.deviance # NULL モデルの deviance
x.glm.nb$df.null # NULL モデルの自由度
pchisq( x.glm.nb2$null.deviance , x.glm.nb$df.null , lower = F)
* residual の種類，data$residuals と residual関数は working 2006 07 31 [#id512ced]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p. 124
response, pearson, working とあるが，
解析結果の $residuals で出力されるのは working
その他の residual を得るには関数 residual を使う．
* 警告メッセージを抑制する 2006 08 01 [#qdb0ef21]
suppressWarnings 関数の引数として，命令を実行する．
* 分布の適合度を一般化線形モデルで行う 2006 08 02 [#q457eec8]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.29 -30.
またカイ自乗分布は，平均 ＝ 自由度，分散 ＝ 2*自由度に従う．
* 飽和モデルの自由度は 0 2006 08 02 [#h3819525]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.33
* predicted と fitted の違い 2006 08 02 [#ba224f21]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.36
* expression の使い方 2006 08 18 [#u90951e9]
描画を行う
si.op1<- sichel.opti$par[1]
si.op2<- sichel.opti$par[2]
text.main<- paste("Sichel's Compound Poisson Distribution\n", text.main)
barplot(both, col = rep(c("white", "green"), length(kekka)), names.arg = cate.label, ylab = "Frequency", xlab = "Cases", main = text.main )
legend(legend.x - 25, legend.y, c("Observed", "Expected"), fill = c("white", "green"))
# legend(legend.x + 15, legend.y - 8, c(expression(paste(alpha, sichel.opti$par[1], beta, sichel.opti$par[2], chi^2, " P (freq > 1) =" ))), chi.score1))
legend(legend.x, legend.y - 8, c(expression(paste(alpha, " = ")), si.op1, expression(paste(beta, " = "))), si.op2, expression(paste(chi^2, " P (freq > 1) =" ))), chi.score1))
* dnorm を使った期待値の計算 2006 08 19 [#t98443bf]
正確に行うには 青木先生のサイト
http://aoki2.si.gunma-u.ac.jp/lecture/GoodnessOfFitness/normaldist.html
http://aoki2.si.gunma-u.ac.jp/R/fit_normal.html
にあるように，各階級の級中心の確率を求め，その階級の期待度数は，上の級の
確率との差にデータ数を乗じることで求まる．
# Crawley, R, p.56
means<- numeric(10000)
for(i in 1:10000){
means[i]<- mean(runif(5) * 10)
}
# sum(means)# 49946.35
length(means)
hist.breaks<- hist(means, ylim=c(0,1600))$breaks
hist.breaks# 1 から 10 の範囲 10個を 2 倍して，0.5, 1, ... としている．
# 従って，length(means) の半分を dnorm の計算結果に乗じる
hist(means, ylim=c(0,1600))
yv lines(xv, yv)
############################### これを形式的に行うには #############
length(hist.breaks)
length(dnorm(hist.breaks, mean = mean(means), sd = sd(means)))
(dnorm.test<- sum(dnorm(hist.breaks, mean = mean(means), sd = sd(means))) * length(means))
sum(means)# mean(means) * length(means)
# 従って
# length(dnorm(hist.breaks, mean = mean(means), sd = sd(means)))
sum(means)#
sum(dnorm.test<- sum(dnorm(hist.breaks, mean = mean(means), sd = sd(means))) * length(means))
# sum(dnorm.test<- sum(dnorm(hist.breaks, mean = mean(means), sd = sd(means))) * sum(means))
test.yv<- dnorm(hist.breaks, mean = mean(means), sd = sd(means)) * length(means) / sum (dnorm(hist.breaks, mean = mean(means), sd = sd(means)) * length(means)) / sum(means)#
test.yv<- dnorm(hist.breaks, mean = mean(means), sd = sd(means)) * length(means) * sum (dnorm(hist.breaks, mean = mean(means), sd = sd(means)) * length(means)) / sum(means)#
hist(means, ylim=c(0,1600))
lines(hist.breaks, test.yv)
##########################################
# 実測値
# log.bread<- seq(min(log.V1) - 0.5, max(log.V1) + 0.5, 0.25)
log.bread<- seq(min(log.V1), max(log.V1)+1, 0.25)
log.V1.ob<- as.vector(table(cut(log.V1, breaks = log.bread)))
length(log.V1.ob)
sum(log.V1.ob)
# 期待値
log.V1.dnorm<- dnorm(log.bread[-(length(log.bread))], mean = mean(log.V1), sd = sd(log.V1))
log.V1.exp<- log.V1.dnorm * length(log.V1) / (sum(log.V1.dnorm * length(log.V1)) / sum(log.V1.ob) )
sum(log.V1.exp)
* merge 関数でデータフレームに要素を追加 2006 09 16 [#y21bc4c8]
test<- data.frame(x = 1, y = 2)
> test
x y
1 1 2
> merge(test, data.frame(x = 3, y = 4), all = T)
x y
1 1 2
2 3 4
* リストの中身を全て取り出す方法 2006 09 16 [#w7ba9652]
################ リストの中身を全て取り出す
x.test<- 1:5
y.test<- 6:8
z.test<- 9:12
xyz<- list(x.test, y.test, z.test)
for(i in 1:length(xyz)){
for(j in 1:length(xyz[[i]])){
print (xyz[[i]][j])
}
}
* if()else{} でelse で始まる文は認められない 2006 09 18 [#ie024e04]
if(bun.p >= 4){
ex.freq<- 5
}else{
ex.freq<- 1
}
は良いが
if(bun.p >= 4){
ex.freq<- 5
}else{
ex.freq<- 1
}
はエラーとなる．
* exists 関数でオブジェクトを確認 2006 09 18 [#m5673588]
if(!exists("myobject"))
* 自作プログラム，観測データとその期待値の区間を調整して，頻度をはかる 2006 09 18 [#x097fc30]
~/daigaku/GakubuKeihi/bun.cut.R より
########## 頻度表の区間をまとめる
# オリジナル
#
bun.df0<- bun.df[-1,]# ゼロ頻度をけずる
#
length(bun.df0$freq)
# bun.p<- 3
########### 指定された区間でデータをまとめる
# bun.p が 1 ならば，元データと同じこと
bun.cut<- data.frame(cate = 1, freq = 0)
xx<- 0
for(i in 1:length(bun.df0$freq)){
xx<- xx + bun.df0[i,2]
z1<- ceiling(i / bun.p)
z2<- i %% bun.p
if(z2 == 0) {
bun.cut[z1,1]<- z1
bun.cut[z1,2]<- xx
xx<- 0
}
if((z2 != 0) && (i == length(bun.df0$freq))){
bun.cut[z1,1]<- z1
bun.cut[z1,2]<- xx
}
}
length(bun.cut$freq)
##########################
# 期待値 の頻度を指定された区間でまとめる
#
theo.cut<- data.frame(cate = 1, freq = 0)
xx<- 0
for(i in 1:length(theo.bun.0)){
xx<- xx + theo.bun.0[i]
z1<- ceiling(i / bun.p)
z2<- i %% bun.p
if(z2 == 0) {
theo.cut[z1,1]<- z1
theo.cut[z1,2]<- xx
xx<- 0
}
if((z2 != 0) && (i == length(bun.df0$freq))){
theo.cut[z1,1]<- z1
theo.cut[z1,2]<- xx
}
}
####################################
# 期待値が ex.freq 未満のセルをまとめる
# ex.freq<- 5
theo.under<- data.frame(cate = 0, freq = 0)
un<- 1
xx<- 0
yy<- numeric(0)
zz<- list(0)
for(i in 1: length(theo.cut$freq)){
if(theo.cut[i,2]< ex.freq){
yy<- c(yy,i)# ex.freq 以下の行番号を記録
# print(yy)
xx<- xx + theo.cut[i,2]
if(xx >= ex.freq){
zz[[un]]<- yy #マージが実行された行を記録
yy<- numeric(0)
theo.under[un,1]<- i #マージ実行時のカテゴリを記録
theo.under[un,2]<- xx ##マージによる融合合計頻度
un<- un + 1
xx<- 0
}
}
else {
if((xx != 0) && xx< ex.freq){# 頻度は5以上だが，前のセルが残っている
zz[[un]]<- c(yy,i) #マージが実行された行を記録
yy<- numeric(0)
theo.under[un,1]<- i #マージ実行時のカテゴリを記録
theo.under[un,2]<- xx + theo.cut[i,2]# #マージによる融合合計頻度
un<- un + 1
xx<- 0
}
}
###### 末尾に達した
if(i == length(theo.cut$freq) ){
if(xx > 0){# そして最後の行のセルが記録されていない
if(un == 1) {# ここまで指定頻度以下のセルが無かった場合
theo.under[un,1]<- i
theo.under[un,2]<- xx
zz[[un]]<- i
}else{# 5 以下の頻度であるが，そのまま出力する#最後は5以下でもよしとする
theo.under[un,1]<- i
theo.under[un,2]<- xx
zz[[un]]<- yy
}
}
}
}
###併合結果をマージする
if(zz[[1]][1] != 0){#併合が行われているならば
zz.vec<- unlist(zz)
theo.target<- theo.cut[-zz.vec,]
bun.target<- bun.cut[-zz.vec,]
if(nrow(theo.target) == 0){#全部が併合対象になっている場合
theo.target<- theo.cut
bun.target<- bun.cut
}
#######
for(i in 1:length(zz)){
theo.target<- merge(theo.target, theo.under[i,], all = T)
bun.target<- merge(bun.target, data.frame(cate = theo.under[i,1],
freq = sum(bun.cut[zz[[i]],2])), all = T)
}
}else{
theo.target<- theo.cut
bun.target<- bun.cut
}
##### 以下三つの出力 (文の数) は等しくなければならない
sum(theo.target$freq)
sum(bun.target$freq)
length(buntyo$V1)
* 0-負の二項分布のパラメータ推定 2006 09 29 [#g3fecf27]
source("/home/ishida/research/statistics/wyshak.R")
あらかじめ次の計算をしておくこと
なお nb0.table ではゼロ頻度のカテゴリも登録しておくこと
(nb0.mean<- mean(nb0.data)) # mean(nb0.table) ではない
(nb0.sum<- sum(nb0.table)) # sum(nb0.data) ではない
(nb0.var<- var(nb0.data))
(nb0.ratio<- nb0.table[1] / sum(nb0.table))#length(nb0.data)
(w0<- nb0.mean / nb0.var * ( 1- nb0.ratio))
(k0<- (w0 * nb0.mean - nb0.ratio) / (1 - w0))
################################ Wyshak
## ################################ Wyshak
## author = {G. Wyshak},
## title = {A Program for Estimating the Parameters of the Truncated Negative Binomial Distribution},
## journal = {Journal of the Royal Statistical Society. Series C, Applied statistics},
## start
## Brass のデータ
## start
nb0.data<- c(rep(1,49), rep(2,56), rep(3,73), rep(4,41), rep(5,43), rep(6,23), rep(7,18), rep(8,18), rep(9,7), rep(10,7), rep(11,3), rep(12,2))
(nb0.mean<- mean(nb0.data)) # mean(nb0.table) ではない
(nb0.sum<- sum(nb0.table)) # sum(nb0.data) ではない
(nb0.var<- var(nb0.data))
(nb0.ratio<- nb0.table[1] / sum(nb0.table))#length(nb0.data)
(w0<- nb0.mean / nb0.var * ( 1- nb0.ratio))
(k0<- (w0 * nb0.mean - nb0.ratio) / (1 - w0))
########################################################## end
## Sampford のデータ
## start
nb0.data<- c(rep(1,11), rep(2,6),rep(3,4), rep(4,5), 6, rep(8,2), 9, 11, 13)
(nb0.table<- table(nb0.data))
# 0 頻度のカテゴリを登録する
#
nb0.df<- data.frame(nb0.table)
y<- 1:max(nb0.data)
nb0.df2<- data.frame(cate = y, freq = c(rep(0, length(y))))
z<- 0
# nb0.df2[1,] = c(0,0)
for(i in 1:nrow(nb0.df2)){
if(nb0.df2[i,]$cate == nb0.df[i-z,]$nb0.data){
nb0.df2[i,]$freq<- nb0.df[i-z,]$Freq
}
else{
z<- z + 1;
next;
}
}
(nb0.mean<- mean(nb0.data)) # mean(nb0.table) ではない
(nb0.sum<- sum(nb0.df2$freq)) # sum(nb0.data) ではない
(nb0.var<- var(nb0.data))
(nb0.ratio<- nb0.df2[1,]$freq / sum(nb0.df2$freq))#length(nb0.data)
(w0<- nb0.mean / nb0.var * ( 1- nb0.ratio))
(k0<- (w0 * nb0.mean - nb0.ratio) / (1 - w0))
nb0.table<- nb0.df2$freq
####################################################### end
## start
wyshak.wk nw<- (nb0.sum *z[2] * log(z[1])) - (nb0.sum * log(1-z[1]^z[2]))
r<- length(nb0.table)# 頻度表 途中のゼロ頻度を省略してはいけない
nx<- 0
for(j in 1:r){
nx<- nx + j * nb0.table[j]# * log(1 - z[1])
}
ny<- 0
for(j in 1:r){
ny<- ny + nb0.table[j] * log(factorial(j))
}
nz<- 0
for(j in 1:r){
nz2<- 0
for(j2 in 1:j){
nz2<- nz2 + log(z[2] + j2 - 1)
}
nz<- nz + nb0.table[j] * nz2
}
return( nw + nx * log(1 - z[1]) - ny + nz )
}
#
optim(c(w0, k0), wyshak.wk, control = list(fnscale = -1))# 最大化する
system.time(optim(c(w0, k0), wyshak.wk, control = list(fnscale = -1)))
* optimize の使いかた 2006 11 08 [#a7d9a62a]
渡部『ベイズ統計学入門』 p.72
モデル分布がベルヌーイ分布である母集団から n = 3 の標本を抽出したとき，その観測値は
c(1,0,1)
であった．母数の最尤推定値を求めよ．
f<- function(x) x^2 * (1-x)
optimize(f,c(0.001, 0.999),maximum=T)
* fitdistr の使いかた 2006 11 08 [#f68b83ea]
http://www.okada.jp.org/RWiki/index.php?R%A4%CE%C5%FD%B7%D7%B2%F2%C0%CF%B4%D8%BF%F4Tips
より
library(MASS)
> mydt<- function(x, m, s, df) dt((x - m)/s, df)/s # 密度関数を独自定義
> fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower = c(-Inf, 0))
m s # 二つのパラメータの推定値と推定標準偏差
-0.01069635 1.04409435
( 0.07222562) ( 0.05434249)
> set.seed(123) # 乱数種
> x4<- rnegbin(500, mu = 5, theta = 4) # 負の二項分布に従う乱数500個発生
> fitdistr(x4, "Negative Binomial") # 負の二項分布を当てはめ
size mu # 二つのパラメータの推定値と推定標準偏差
4.2178713 4.9439803
(0.5048242) (0.1465538)
Warning messages: # 関連警告
* poisson 分布で予測する 2006 11 08 [#b9f98144]
x<- c(0,0,0,0,0,2)
library(MASS)
fitdistr(x, "Poisson")
## lambda
## 0.3333333
## (0.2357023)
# 0 人現れる確率
exp(-0.3333) * 0.3333^0 / factorial(0)
#[1] 0.7165552
# 1 人現れる確率
0.3333/(0+1) * 0.7165552
# [1] 0.2388278
# 2 人現れる確率
0.3333/(1+2) * 0.2388278
* 多元分割表への対数線形モデル 2006 11 23 [#f7234641]
&aname(matsuda121);
質的情報の多変量解析 (松田 紀之 著)# p.121
matsuda.121<- array(c(32,86, 11, 35, 61,73,41,70),
dim = c(2,2,2),
dimnames = list(height=c("over","lower"),
diam = c("lower","heigher"),
species = c("sagrei", "distichus")))
class( matsuda.121 )<- "table"
matsuda.121
# p.129
matsuda.121.loglm<- loglm(~ height * diam * species, data = matsuda.121 )
stepAIC(matsuda.121.loglm )
Step: AIC= 14.03
~height + diam + species + height:species + diam:species
Df AIC
<none> 14.026
- height:species 1 22.431
- diam:species 1 24.632
Call:
loglm(formula = ~height + diam + species + height:species + diam:species,
data = matsuda.121, evaluate = FALSE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 2.025647 2 0.3631921
Pearson 2.017364 2 0.3646994
同じことだが
matsuda.121.loglin<- loglin( matsuda.121, list(c(1,3), c(2,3)) )
matsuda.121.loglin
# ↑ c(1,2) を外してよいか．つまり heigh と diam は独立だとみなし，無視してよいか
# 実際に外した結果を，飽和モデルと比較する
# 0.3631921 なので，外してよい ↓
1 - pchisq(matsuda.121.loglin$lrt, matsuda.121.loglin$df)
1 - pchisq(matsuda.121.loglin$pearson, matsuda.121.loglin$df)#
よって species さえ分かれば，height や diam の間に特に関係は無い．
* ホームディレクトリへのパッケージのインストール 2006 12 02 [#jc18d01d]
次のように，デフォルト以外のライブラリの場所を確認
Sys.getenv("R_LIBS")
必要があれば
Sys.setenv(R_LIBS="C:/R")
install.packages("MyPackages", lib = "C:/R") # 第二引数はオプション
library("MyPackages")
* 日本語フォントの利用 2006 12 05 [#s46937ef]
http://www.okada.jp.org/RWiki/index.php?%C6%FC%CB%DC%B8%EC%B2%BD%B7%C7%BC%A8%C8%C4 より
フォントの確認方法
postscriptFonts()
postscriptFonts()$Japan1
options(X11fonts = c("-alias-gothic-%s-%s-*-*-%d-*-*-*-*-*-*-*",
"-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
あるいは
options(X11fonts = c("-alias-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*","-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
ps.options(family= "Japan1GothicBBB")
あるいは
ps.options(family= "Japan1Ryumin")
Rを起動して, 上の2行をカットアンドペーストして
plot(sin,-pi,pi,main="正弦")
とかすると多分,日本語が出てくると思います.
ディスプレイ(x11,pdf,jpeg)のフォントについては help(X11)とかが参考になるでしょうし,
PDFは help(pdf) を見てください.
デフォルトとしてフォントを設定するなら help(Rprofile) を見てください. -- なかま 2006-11-24 (金) 22:14:57
あるいは
日本語表示をするために、 X11デバイスへの出力には
> options(X11fonts = c("-*-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*",
"-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
options(X11fonts = c("-alias-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*",
"-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
としておいて、
> plot(適当な引数)
日本語表示のためにpostscript(EPS画像ファイル)への出力には、
> postscript("output.eps",family="Japan1Ryumin")
としておいて、
> plot(適当な引数)
> dev.off()
* Fonts の設定 X11fonts [#a7f9c156]
R のデフォルトは
> options()$X11fonts
[1] "-adobe-helvetica-%s-%s-*-*-%d-*-*-*-*-*-*-*"
[2] "-adobe-symbol-medium-r-*-*-%d-*-*-*-*-*-*-*"
Fedora core5 では以下のXLFD 情報などがきれい．
options(X11fonts = c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*",
"-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*-0"))
options(X11fonts = c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*-*",
"-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
options(X11fonts = c("-shinonome-mincho-%s-%s-*-*-%d-*-*-*-*-*-*-*",
"-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*"))
ただし，座標ラベルが指定されていない場合，無意味な文字が出力される．
x<- c("あ","い","う","え","お")
y<- c("ア","イ","ウ","エ","オ")
plot(1:5, 1:5, xlab = "これはテスト", type = "n", ylab = "")
text(1:5, 1:5, x, cex = 0.7)
legend("topright" , legend = c("日本語ひらがら"))
RWiki には以下の情報があった．
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=R%B7%C7%BC%A8%C8%C4&word=X11fonts
X11のフォントの設定は.Rprofileに以下を加えて下さい
options(X11fonts=c("-kochi-mincho-%s-%s-normal--%d-*-*-*-*-*-*-*",
"-adobe-symbol-medium-r-*-*-%d-*-*-*-*-*-*-*" ))
X11日本語フォントの一覧は,
xlsfonts -fn "*jis*"
で得られます.
全プラットフォーム共通
日本語のPostScript?及びPDFを扱う場合は, .Rprofileに以下を記して下さい
setHook(packageEvent("grDevices", "onLoad"),
function(...) grDevices::ps.options(family="Japan1"))
PostScript?,PDFのフォントの一覧を得るならコマンドラインより
postscript
cbind(lapply(postscriptFonts(),function(x){x$family}))
pdf
cbind(lapply(pdfFonts(),function(x){x$family}))
尚,イラストレータ等外部にデータを渡す場合は日本語フォントの場合はEPSでは無く
PDFをお薦めします
-- なかま 2006-04-25 (火) 00:40:39
* データフレームから，列名を指定して抽出 2006 12 07 [#xf85f0b8]
例えば[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.67 で
data(skulls, package = "HSAUR")
means<- aggregate(skulls[, c("mb", "bh", "bl", "nh")],
list(epoch = skulls$epoch), mean)
# 他に
cov(iris[, -5]) # 共分散行列．var(iris[,-5] でも同じ
var(iris[, colnames(iris) != "Species"]) # 共分散行列．条件指定を変更
cor(iris[sapply(iris, is.numeric)]) # 相関行列 # 条件指定を変更
* Contrasts Matrix の意味 2006 12 08 [#l8cd5384]
&aname(ContrastsMatrix);
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.75 下
* 外れ値を表示する方法 2006 12 08 [#jd2e37de]
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.78
data("clouds", package = "HSAUR")
layout(matrix(1:2, nrow = 2))
bxpseeding xlab = "Seeding")
bxpsecho xlab = "Echo motion")
rownames(clouds)[clouds$rainfall %in% c(bxpseeding$out, bxpsecho$out)]
* subset の使いかた 2006 12 08 [#lf446f10]
* legend の使いかた 2006 12 08 [#f2a50ee5]
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.51
cmh_test( classification ~ treatment, data = Lanza, scores = list(classification = c(0, 1, 6, 17, 30)),
subset = Lanza$study == "I")
cmh_test( classification ~ treatment, data = Lanza, scores = list(classification = c(0, 1, 6, 17, 30)),
subset = Lanza$study == "II")
# p.82
psymb<- as.numeric(clouds$seeding)
plot(rainfall ~ cloudcover , data = clouds, pch = psymb)
abline(lm(rainfall ~ cloudcover, data = clouds,
subset = seeding == "no"))
abline(lm(rainfall ~ cloudcover, data = clouds,
subset = seeding == "yes"), lty = 2)
legend("topright" , legend = c("No seeding", "Seeding"),
pch = 1:2, lty = 1:2, bty = "n")
* quasilikelihood による overdispersion の考慮 2006 12 13 [#o2919a11]
Everit & Hothorn p.106
また overdispersion が生じるのは，しばしばデータが独立でないため
* mixture model 2006 12 15 [#bb74ed7e]
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p.117
* grep 関数の使いかた 2006 12 22 [#od716620]
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] p. 166
例えば以下のコードは
BtheB[, grep("bdi", names(BtheB))]
BtheB データフレームで，列名に "bdi" がついている列番号を抽出している．
さらに以下は
subset(BtheB, treatment == "TAU")[, grep("bdi", names(BtheB))]
BtheB データで，列 treatment が "TAU" であるサブセットから，
"bdi" を名前として含む列を取り出している．
* reshape 関数の使いかた 2006 12 22 [#m7ddbbce]
[[Everitt & Hothron, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] chapter 10 p. 165
data("BtheB", package = "HSAUR")
BtheB$subject<- factor(rownames(BtheB))
nobs<- nrow(BtheB)
BtheB.long<- reshape(BtheB, idvar = "subject",
varying = c("bdi.2m", "bdi.4m", "bdi.6m", "bdi.8m"), direction = "long")
# 石村貞夫「分散分析のはなし」の例．データは p.77
# データフレームの作成
# kaeru<- read.csv("africanFlog.csv",header=F)
# 始めに横長形式のデータとして与えられていたとする．
frogs<- data.frame(stage = c("A1", "A2", "A3", "A4", "A5"),
sample1 = c(12.2, 22.2, 20.8, 26.4, 24.5),
sample2 = c(18.8, 20.5, 19.5, 32.6, 21.2),
sample3 = c(18.2, 14.6, 26.3, 31.3, 22.4))
frog.res<- reshape(frogs, idvar = "stage",
varying = list(names(frogs[2:4])),v.names = "data", direction = "long")
# ここで作成される time 変数は，各数値の順番を表す
# (時間ごとにデータを取っている場合などに重要)
frog2<- frog.res[order(frog.res$stage),] # 並び替え
frog2
* 判別分析での適合率 2006 12 26 [#n9eb6c77]
http://biking.taiiku.tsukuba.ac.jp/wiki/index.php?%C5%FD%B7%D7%A4%A2%A4%EC%A4%B3%A4%EC
判別分析では、モデルに含まれるパラメーター数が多くなればなるほど、
あてはめ誤差（残差平方和）は小さくなるということがあります。
そこで、単に残差平方和の大小を比較するだけでなく、パラメータ数も考慮したモデルの選択を行います。
その際の基準となるものがAICです。wleライブラリ中のmle.aic関数を使ってモデルを選択します。
>library(MASS)
>library(wle)
>result<-mle.aic(lda(class~.,data=DATA))
>summary(result,num=20)
mle.cvでやった場合にも判別率は出力されませんでした。AIC値の代わりにcv値が出力されていますが、どういう値なのかよくわかりません。
上記に続いて、判別関数を求める手順、交差妥当性（クロス・バリデーション）のチェックをやって判別率を比べる手順を示します。
> (z<-lda(class~.,data=DATA))
>apply(z$means%*%Z$scaling,2,mean)
この結果、それぞれの変数の係数と定数が出力されます。
>predict(z)$x
>predict(z)$class
>predict(z)$posterior
と入力すればこの判別関数を用いた判別得点、各個体が判別されたグループのラベル、
各個体がどのグループに判別されているかに関する事後確率（0から1）が出力されます。
判別率を計算するためのクロス表の出力は次のようにします。
>table(data[,],predict(z)$class)
交差妥当性のチェックを行ったクロス表は
>DATA.CV<-lda(class~.,data=DATA,CV=T)
>(lda.tab<-table(DATA[,],data=DATA.CV$class))
で出力できます。判別率、誤判別率はそれぞれ
>sum(lda.tab[row(lda.tab)==col(lda.tab)])/sum(lda.tab)
>sum(lda.tab[row(lda.tab)!=col(lda.tab)])/sum(lda.tab)
交差確認（cross validation）はデータを2分し,一方のサンプルで判別式をつくりその判別式で残りのサンプルを判別して,
判別式の有効性を見る方法です。データを2分する以外に3,4,5,10分する方法もありますが、
どれを用いるかは明確な基準がなくデータサイズに依存するそうです。
この特殊なケースとして1 つ取っておき法（reave-one-out cross-validation）があります。
データから1つの個体を取り除いて判別分析を行い、取り除かれた個体で判別モデルの評価を行うという作業を全ての個体に繰り返して行います。
「R」でこれが出来るということなので紹介します。
ここでは「R」の中にあるirisデータを例に使います。
> library(MASS)
> iris.CV<-lda(Species~.,data=iris,CV=T)
> (lda.tab<-table(iris[,5],iris.CV$class))
デフォルトではCV=FALSEとなっているのでCV=Tにすると1つ取っておき法（reave-one- out cross-validation）ができます。
ちなみに判別関数lda(Linear Discriminant Analysis),qda(Quadratic Discriminant Analysis)には1つ取っておき法（reave-one-out cross-validation）以外の交差確認機能は用意されていないそうです。
詳しい説明は,同志社大学Jin先生のページhttp://www1.doshisha.ac.jp/~mjin/R/の「Rと判別分析」を見ると書いてあります。
* １つ取っておき法とjackknife法はほとんど同義語として使われるようです。
しかし，jackknife法よりも良さそうな方法，bootstrap （ブートストラップ）法があります。
データからm個の標本を復元抽出（sampling wtih replacement）し，それらの標本から繰り返しパラメータ推定値を求める方法です。
これを線型モデル（回帰分析）で使った例が，ヴェナブルズとリプリーの「S-PLUSによる統計解析」シュプリンガー・フェアラーク東京(2001)の211〜214ぺージにあります。
* bootstrap法をつかって判別分析の交差確認を実行するには，
Rのipredライブラリのerrorest()関数が使えます。ipredはimproved predictorsの略。
> library(ipred)
> mypredict.lda<- function(object, newdata) predict(object, newdata=newdata)$class
> errorest(class ~ ., data=soccer, model=lda, estimator="boot", predict=mypredict.lda)
mypredict.ldaはclass labelだけを返させるのに必要な手続きらしい。
* 判別分析での変数選択 2006 12 26 [#qa081b4c]
- 線型判別分析
klaRライブラリ中のstepclass関数を使う。
methodに線型判別分析関数のldaを与える。ステップワイズの方向は変数増減がデフォルトとなっている。
library(klaR); stepclass(Y ~ ., data=DATA, method="lda")
- ロジスティック回帰
MASSライブラリ中のstepAIC関数を使ってステップワイズ変数選択（変数増減法がデフォルト）をおこなう。
引数には一般化線型モデルの関数glmの応答分布の関数に二項分布（binomial）の結果を与える（logitがリンク関数のデフォルトとなっている）。
library(MASS); stepAIC(glm(Y ~ ., data=DATA, family=binomial))
* データフレームと文字列 2006 12 27 [#z47228dd]
test<- data.frame(v1 = c("A","B","C"), v2 = c(1,2,3))
test$v1[test$v1 == test$v1[1]]<- "AA"
これは
Warning message:
無効な因子水準です。NA が発生しました in: `[<-.factor`(`*tmp*`, iseq, value = "AA")
となる．そこで，
test$v1<- as.character(test$v1)
is.factor(test$v1); FALSE
test$v1[test$v1 == test$v1[1]]<- "AA"
とすれば，交換はできる
test<- data.frame(v1 = c("X","Y","Z","X","Y"), v2 = c(1,2,3,4,5))
test$v1<- as.character(test$v1)
test.lab<- c("X","Y","Z")
for(i in 1:length(test.lab)){
test$v1[test$v1 == test.lab[i]]<- LETTERS[i]
}
* parse eval の使いかた 2006 12 27 [#tcf54943]
もしも列名などで，ループ処理したいときは，以下のようにする．
test<- data.frame(v1 = c("A","B","C"), v2 = c(1,2,3))
sum( eval (parse(text = paste("test$v", 2, sep="") ) ))
[1] 6
ちなみに，以下は動作しない
xxx<- "v2"
sum( eval (parse(text = "test$xxx") ))
ただし，一般には
sum(text[,i])
として列を指定してループ処理する．
他に /daigaku/kakei/forDrKakei.R を参照のこと
* mtext の使いかた 2006 12 28 [#r9850e3f]
「Rの基礎とプログラミング技法」（シュプリンガージャパン）169ページのグラフィックス
~/research/statistics/ligges.R
par(oma=rep(3, 4), bg="white")
plot(c(0, 1), c(0, 1), type="n", ann=FALSE, axes=FALSE)
rect(-1, -1, 2, 2, col="white")
box("figure")
par(xpd=FALSE)
rect(-1, -1, 2, 2, col="white")
box("plot", lty="dashed")
text(.5, .5, "plot region", cex = 1.6)
mtext("figure region", side=3, line=2, adj = 1, cex = 1.4)
for (i in 1:4)
mtext(paste("inner margin", i), side=i, line=1, outer= FALSE)
for (i in 1:4)
mtext(paste("outer margin", i), side=i, line=1, outer=TRUE)
mtext("device region", side=3, line=2, outer=TRUE, adj = 1, cex = 1.2)
box("outer", col="black")