Site admin: motohiro PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 7.4.33
HTML convert time: 0.148 sec.
個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-m.at.ias.tokushima-u.ac.jp. までご連絡ください． 個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-m(この部分を"@"に変更下さい)ias.tokushima-u.ac.jp までご連絡ください． [[R_fromOldHtml2]] #contents * バートレットの多群の等分散性 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))) ) 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]]を参考にさせていただきました．以下引用です！(石田)． 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) 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 [[このサイト: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 ''以下は 次からの引用です！ 石田'' ''以下は[[このサイト: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) 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) 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) 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 Wood, GLM, p.33 モデルに不必要な説明変数が多く存在する * A % in % B による条件判定とデータ抽出 2006 04 28 [#ucc36199] Wood, p.35 Wood, GLM, p.35 sperm.comp1$m.vol<- 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 他に [[Everitt & Hothorn, 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 Wood, GLM, 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 上 [[Everitt & Hothorn, 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 Everitt & 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 Everitt & 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 村上征勝先生の「工業統計学」, 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.）。 * 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 下 # 間瀬先生著 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. [[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-Generalized-Nonparametric/dp/158488424X/]] p.120 (さらに 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) # 石田 追加，これは便宜上の説明変数を作ったに過ぎない 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) 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) 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)) 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 Crawley, Statistics using 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, ... としている． 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))) 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), 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), # 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, 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) - 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) 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以下でもよしとする # 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") source("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)) 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 ############### 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) f <- function(x) x^2 * (1-x) optimize(f,c(0.001, 0.999),maximum=T) * fitdistr の使いかた 2006 11 08 [#f68b83ea] [[RjpWikiを参照: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 mydt<- function(x, m, s, df) dt((x - m)/s, df)/s # 密度関数を独自定義 > fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, fitdistr(x2, mydt, list(m = 0, s = 1), df = 9, lower = c(-Inf, 0)) m s # 2パラメータの推定値と推定標準偏差 -0.01069635 1.04409435 ( 0.07222562) ( 0.05434249) > set.seed(123) # 乱数種 > x4<- rnegbin(500, mu = 5, theta = 4) set.seed(123) # 乱数種 x4<- rnegbin(500, mu = 5, theta = 4) # 負の二項分布に従う乱数500個発生 > fitdistr(x4, "Negative Binomial") # 負の二項分布を当てはめ size mu # 2パラメータの推定値と推定標準偏差 4.2178713 4.9439803 (0.5048242) (0.1465538) fitdistr(x4, "Negative Binomial") # 負の二項分布を当てはめ size mu # 2パラメータの推定値と推定標準偏差 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 人現れる確率 # 1 人現れる確率 0.3333/(0+1) * 0.7165552 # [1] 0.2388278 # 2 人現れる確率 0.3333/(1+2) * 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 # 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] [[RjpWikiを参照: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 * Fonts の設定 X11fonts [#a7f9c156] R のデフォルトは > options()$X11fonts [1] "-adobe-helvetica-%s-%s-*-*-%d-*-*-*-*-*-*-*" [2] "-adobe-symbol-medium-r-*-*-%d-*-*-*-*-*-*-*" 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-*-*-*-*-*-*-*", options(X11fonts = c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*", "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*-0")) options(X11fonts = c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*-*", c("-shinonome-gothic-%s-%s-normal--%d-*-*-*-*-*-*-*-*", "-adobe-symbol-*-*-*-*-%d-*-*-*-*-*-*-*")) options(X11fonts = c("-shinonome-mincho-%s-%s-*-*-%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]] * データフレームから，列名を指定して抽出 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 で 例えば[[Everitt & Hothorn, 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 下 [[Everitt & Hothorn, 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 [[Everitt & Hothorn, 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)] 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 [[Everitt & Hothorn, 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") 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") 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 [[Everitt & Hothorn, A Handbook of Statistical Analyses Using R:http://www.amazon.co.jp/Handbook-Statistical-Analyses-Using-R/dp/1584885394/]] 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 [[Everitt & Hothorn, 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 [[Everitt & Hothorn, 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 [[Everitt & Hothorn, 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)) 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 変数は，各数値の順番を表す # (時間ごとにデータを取っている場合などに重要) # ここで作成される time 変数は，各数値の順番を表す # (時間ごとにデータを取っている場合などに重要) frog2<- frog.res[order(frog.res$stage),] # 並び替え frog2 * 判別分析での適合率 2006 12 26 [#n9eb6c77] 以下は [[Rjp Wiki:http://biking.taiiku.tsukuba.ac.jp/wiki/index.php?%C5%FD%B7%D7%A4%A2%A4%EC%A4%B3%A4%EC]] より''の引用です！ 石田'' 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) 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) (z<-lda(class~.,data=DATA)) apply(z$means%*%Z$scaling,2,mean) この結果、それぞれの変数の係数と定数が出力されます。 >predict(z)$x >predict(z)$class >predict(z)$posterior predict(z)$x predict(z)$class predict(z)$posterior と入力すればこの判別関数を用いた判別得点、各個体が判別されたグループのラベル、各個体がどのグループに判別されているかに関する事後確率（0から1）が出力されます。判別率を計算するためのクロス表の出力は次のようにします。 >table(data[,],predict(z)$class) table(data[,],predict(z)$class) 交差妥当性のチェックを行ったクロス表は >DATA.CV<-lda(class~.,data=DATA,CV=T) >(lda.tab<-table(DATA[,],data=DATA.CV$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) 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)) > 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と判別分析」を見ると書いてあります。 1つ取っておき法と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) 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関数を使う。mehodに線型判別分析関数の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] 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 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") 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")