トップ
新規
一覧
単語検索
最終更新
ヘルプ
ログイン
アールメカブ
R_fromOldHtml2_2
をテンプレートにして作成
開始行:
個人メモを整理中です．不適切な記述が多々あるかと思います...
[[R_fromOldHtml2]]
#contents
* バートレットの多群の等分散性 bartlett.test 2006 04 01 [...
永田靖著『統計的多重比較法の基礎』サイエンティスト社
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/basic...
対数オッズは，通常「ロジット 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-...
* アルファベット順以外のカテゴリ分け関数 gl 2006 04 05 [...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
カテゴリの基準をアルファベット順以外で構成する時などに使う
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-...
> 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 が基準となっ...
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-...
> 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 の別の...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
> 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-...
関数 svd の対象は residuals
* plot 関数への asp 引数 2006 04 05 [#cef1e350]
[[このサイト:http://www.ma.hw.ac.uk/ams/Rhelp/library/gra...
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-...
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-...
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.va...
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-...
cutinc さらに ここ も参照
* log による summary 出力を antilog (exp) 化して解釈する...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
> 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...
帰無仮説どおりのランダムサンプリングによって得られたデー...
ところが、どんな検定でも検定統計量のP-valueを計算すると、...
つまり、 「P-value とは、その値が小さければ小さいほど仮説...
どんな帰無仮説に基いてどんな検定統計量を用いるどんな検定...
* Wald test と，その問題点 2006 04 08 [#m78ef4d2]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* データを中心化 (center) する 2006 04 15 [#f835bdd4]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* iid 独立同型分布 2006 04 15 [#z3694c79]
iid は独立同型分布(independently and identically distribu...
* expand.grid 変数の値を変化させる(ただし他変数は固定) 2...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
expand.grid(temp = seq(-3, 3, 0.1), ibh = 0, ibt = 0)
ibh, ibt の値は平均(標準化された平均値)に固定し，temp を[...
* 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 [#b3ea...
Wood, GLM, p.33
モデルに不必要な説明変数が多く存在する
* A % in % B による条件判定とデータ抽出 2006 04 28 [#uc...
Wood, GLM, p.35
sperm.comp1$m.vol <-
sperm.comp2$m.vol
[sperm.comp2$pair %in% sperm.comp1$subject]
他に [[Everitt & Hothorn, A Handbook of Statistical Analy...
offset の意味 2006 04 28
Wood, GLM, p.45
モデルの含める際，そのパラメータが 1 と仮定する。
* layout 関数によるプロット座標の分割 2006 05 04 [#hb41...
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 & Hothorn, A Handbook of Statistical Analyses U...
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
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/mj...
* bubble plot の効果的な使い方 2006 05 12 [#uc43378a]
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 を抽出 2...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
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 の調整 20...
# 高頻度では、データ区間が飛び飛びになっているので、こ...
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 comm...
* 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.76...
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/...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
このような記事もあった
また次のような記事もあった
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@c...
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...
This result is borne out by simulations. Is this corr...
# 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関数は wor...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
response, pearson, working とあるが，
解析結果の $residuals で出力されるのは working
その他の residual を得るには関数 residual を使う．
* 警告メッセージを抑制する 2006 08 01 [#qdb0ef21]
suppressWarnings 関数の引数として，命令を実行する．
* 分布の適合度を一般化線形モデルで行う 2006 08 02 [#q457...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
またカイ自乗分布は，平均 ＝ 自由度，分散 ＝ 2*自由度に従...
* 飽和モデルの自由度は 0 2006 08 02 [#h3819525]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* predicted と fitted の違い 2006 08 02 [#ba224f21]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* expression の使い方 2006 08 18 [#u90951e9]
描画を行う
si.op1<- sichel.opti$par[1]
si.op2<- sichel.opti$par[2]
text.main<- paste("Sichel's Compound Poisson Distributi...
barplot(both, col = rep(c("white", "green"), length(kek...
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/n...
http://aoki2.si.gunma-u.ac.jp/R/fit_normal.html
にあるように，各階級の級中心の確率を求め，その階級の期待...
確率との差にデータ数を乗じることで求まる．
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, ... として...
# 従って，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.o...
sum(log.V1.exp)
* merge 関数でデータフレームに要素を追加 2006 09 16 [#y...
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...
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"))
* 自作プログラム，観測データとその期待値の区間を調整し...
~/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,], al...
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("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....
(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...
journal = {Journal of the Royal Statistical Society. Se...
## 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...
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]
[[RjpWikiを参照:http:www.okada.jp.org/RWiki/index.php?R%A...
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 # 2パラメータの推定値と推定...
-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 # 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 人現れる確率
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:speci...
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 1...
次のように，デフォルト以外のライブラリの場所を確認
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...
フォントの確認方法
postscriptFonts()
postscriptFonts()$Japan1
* 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...
* データフレームから，列名を指定して抽出 2006 12 07 [#x...
例えば[[Everitt & Hothorn, A Handbook of Statistical Anal...
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 & Hothorn, A Handbook of Statistical Analyses U...
* 外れ値を表示する方法 2006 12 08 [#jd2e37de]
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
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 & Hothorn, A Handbook of Statistical Analyses U...
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 1...
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
また overdispersion が生じるのは，しばしばデータが独立で...
* mixture model 2006 12 15 [#bb74ed7e]
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
* grep 関数の使いかた 2006 12 22 [#od716620]
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
例えば以下のコードは
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 & Hothorn, A Handbook of Statistical Analyses U...
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", "A...
sample1 = c(12.2, 22.2, 20.8, 26.4, 24...
sample2 = c(18.8, 20.5, 19.5, 32.6, 21...
sample3 = c(18.2, 14.6, 26.3, 31.3, 22...
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]
以下は
[[Rjp Wiki:http://biking.taiiku.tsukuba.ac.jp/wiki/index....
判別分析では、モデルに含まれるパラメーター数が多くなれば...
library(MASS)
library(wle)
result<-mle.aic(lda(class~.,data=DATA))
summary(result,num=20)
mle.cvでやった場合にも判別率は出力されませんでした。AIC値...
上記に続いて、判別関数を求める手順、交差妥当性（クロス・...
(z<-lda(class~.,data=DATA))
apply(z$means%*%Z$scaling,2,mean)
この結果、それぞれの変数の係数と定数が出力されます。
predict(z)$x
predict(z)$class
predict(z)$posterior
と入力すればこの判別関数を用いた判別得点、各個体が判別さ...
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,...
この特殊なケースとして1 つ取っておき法（reave-one-out cro...
データから1つの個体を取り除いて判別分析を行い、取り除かれ...
「R」でこれが出来るということなので紹介します。ここでは「...
library(MASS)
iris.CV<-lda(Species~.,data=iris,CV=T)
(lda.tab<-table(iris[,5],iris.CV$class))
デフォルトではCV=FALSEとなっているのでCV=Tにすると1つ取っ...
ちなみに判別関数lda(Linear Discriminant Analysis),qda(Qua...
詳しい説明は,同志社大学Jin先生のページhttp://www1.doshish...
1つ取っておき法とjackknife法はほとんど同義語として使われ...
しかし，jackknife法よりも良さそうな方法，bootstrap （ブー...
データからm個の標本を復元抽出（sampling wtih replacement...
これを線型モデル（回帰分析）で使った例が，ヴェナブルズと...
bootstrap法をつかって判別分析の交差確認を実行するには，
Rのipredライブラリのerrorest()関数が使えます。ipredはimpr...
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に線型判別...
library(klaR); stepclass(Y ~ ., data=DATA, method="lda")
- ロジスティック回帰
-- MASSライブラリ中のstepAIC関数を使ってステップワ...
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])
として列を指定してループ処理する．
* mtext の使いかた 2006 12 28 [#r9850e3f]
「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")
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...
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")
終了行:
個人メモを整理中です．不適切な記述が多々あるかと思います...
[[R_fromOldHtml2]]
#contents
* バートレットの多群の等分散性 bartlett.test 2006 04 01 [...
永田靖著『統計的多重比較法の基礎』サイエンティスト社
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/basic...
対数オッズは，通常「ロジット 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-...
* アルファベット順以外のカテゴリ分け関数 gl 2006 04 05 [...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
カテゴリの基準をアルファベット順以外で構成する時などに使う
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-...
> 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 が基準となっ...
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-...
> 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 の別の...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
> 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-...
関数 svd の対象は residuals
* plot 関数への asp 引数 2006 04 05 [#cef1e350]
[[このサイト:http://www.ma.hw.ac.uk/ams/Rhelp/library/gra...
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-...
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-...
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.va...
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-...
cutinc さらに ここ も参照
* log による summary 出力を antilog (exp) 化して解釈する...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
> 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...
帰無仮説どおりのランダムサンプリングによって得られたデー...
ところが、どんな検定でも検定統計量のP-valueを計算すると、...
つまり、 「P-value とは、その値が小さければ小さいほど仮説...
どんな帰無仮説に基いてどんな検定統計量を用いるどんな検定...
* Wald test と，その問題点 2006 04 08 [#m78ef4d2]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* データを中心化 (center) する 2006 04 15 [#f835bdd4]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* iid 独立同型分布 2006 04 15 [#z3694c79]
iid は独立同型分布(independently and identically distribu...
* expand.grid 変数の値を変化させる(ただし他変数は固定) 2...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
expand.grid(temp = seq(-3, 3, 0.1), ibh = 0, ibt = 0)
ibh, ibt の値は平均(標準化された平均値)に固定し，temp を[...
* 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 [#b3ea...
Wood, GLM, p.33
モデルに不必要な説明変数が多く存在する
* A % in % B による条件判定とデータ抽出 2006 04 28 [#uc...
Wood, GLM, p.35
sperm.comp1$m.vol <-
sperm.comp2$m.vol
[sperm.comp2$pair %in% sperm.comp1$subject]
他に [[Everitt & Hothorn, A Handbook of Statistical Analy...
offset の意味 2006 04 28
Wood, GLM, p.45
モデルの含める際，そのパラメータが 1 と仮定する。
* layout 関数によるプロット座標の分割 2006 05 04 [#hb41...
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 & Hothorn, A Handbook of Statistical Analyses U...
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
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/mj...
* bubble plot の効果的な使い方 2006 05 12 [#uc43378a]
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 を抽出 2...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
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 の調整 20...
# 高頻度では、データ区間が飛び飛びになっているので、こ...
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 comm...
* 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.76...
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/...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
このような記事もあった
また次のような記事もあった
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@c...
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...
This result is borne out by simulations. Is this corr...
# 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関数は wor...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
response, pearson, working とあるが，
解析結果の $residuals で出力されるのは working
その他の residual を得るには関数 residual を使う．
* 警告メッセージを抑制する 2006 08 01 [#qdb0ef21]
suppressWarnings 関数の引数として，命令を実行する．
* 分布の適合度を一般化線形モデルで行う 2006 08 02 [#q457...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
またカイ自乗分布は，平均 ＝ 自由度，分散 ＝ 2*自由度に従...
* 飽和モデルの自由度は 0 2006 08 02 [#h3819525]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* predicted と fitted の違い 2006 08 02 [#ba224f21]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* expression の使い方 2006 08 18 [#u90951e9]
描画を行う
si.op1<- sichel.opti$par[1]
si.op2<- sichel.opti$par[2]
text.main<- paste("Sichel's Compound Poisson Distributi...
barplot(both, col = rep(c("white", "green"), length(kek...
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/n...
http://aoki2.si.gunma-u.ac.jp/R/fit_normal.html
にあるように，各階級の級中心の確率を求め，その階級の期待...
確率との差にデータ数を乗じることで求まる．
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, ... として...
# 従って，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.o...
sum(log.V1.exp)
* merge 関数でデータフレームに要素を追加 2006 09 16 [#y...
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...
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"))
* 自作プログラム，観測データとその期待値の区間を調整し...
~/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,], al...
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("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....
(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...
journal = {Journal of the Royal Statistical Society. Se...
## 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...
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]
[[RjpWikiを参照:http:www.okada.jp.org/RWiki/index.php?R%A...
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 # 2パラメータの推定値と推定...
-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 # 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 人現れる確率
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:speci...
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 1...
次のように，デフォルト以外のライブラリの場所を確認
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...
フォントの確認方法
postscriptFonts()
postscriptFonts()$Japan1
* 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...
* データフレームから，列名を指定して抽出 2006 12 07 [#x...
例えば[[Everitt & Hothorn, A Handbook of Statistical Anal...
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 & Hothorn, A Handbook of Statistical Analyses U...
* 外れ値を表示する方法 2006 12 08 [#jd2e37de]
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
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 & Hothorn, A Handbook of Statistical Analyses U...
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 1...
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
また overdispersion が生じるのは，しばしばデータが独立で...
* mixture model 2006 12 15 [#bb74ed7e]
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
* grep 関数の使いかた 2006 12 22 [#od716620]
[[Everitt & Hothorn, A Handbook of Statistical Analyses U...
例えば以下のコードは
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 & Hothorn, A Handbook of Statistical Analyses U...
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", "A...
sample1 = c(12.2, 22.2, 20.8, 26.4, 24...
sample2 = c(18.8, 20.5, 19.5, 32.6, 21...
sample3 = c(18.2, 14.6, 26.3, 31.3, 22...
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]
以下は
[[Rjp Wiki:http://biking.taiiku.tsukuba.ac.jp/wiki/index....
判別分析では、モデルに含まれるパラメーター数が多くなれば...
library(MASS)
library(wle)
result<-mle.aic(lda(class~.,data=DATA))
summary(result,num=20)
mle.cvでやった場合にも判別率は出力されませんでした。AIC値...
上記に続いて、判別関数を求める手順、交差妥当性（クロス・...
(z<-lda(class~.,data=DATA))
apply(z$means%*%Z$scaling,2,mean)
この結果、それぞれの変数の係数と定数が出力されます。
predict(z)$x
predict(z)$class
predict(z)$posterior
と入力すればこの判別関数を用いた判別得点、各個体が判別さ...
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,...
この特殊なケースとして1 つ取っておき法（reave-one-out cro...
データから1つの個体を取り除いて判別分析を行い、取り除かれ...
「R」でこれが出来るということなので紹介します。ここでは「...
library(MASS)
iris.CV<-lda(Species~.,data=iris,CV=T)
(lda.tab<-table(iris[,5],iris.CV$class))
デフォルトではCV=FALSEとなっているのでCV=Tにすると1つ取っ...
ちなみに判別関数lda(Linear Discriminant Analysis),qda(Qua...
詳しい説明は,同志社大学Jin先生のページhttp://www1.doshish...
1つ取っておき法とjackknife法はほとんど同義語として使われ...
しかし，jackknife法よりも良さそうな方法，bootstrap （ブー...
データからm個の標本を復元抽出（sampling wtih replacement...
これを線型モデル（回帰分析）で使った例が，ヴェナブルズと...
bootstrap法をつかって判別分析の交差確認を実行するには，
Rのipredライブラリのerrorest()関数が使えます。ipredはimpr...
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に線型判別...
library(klaR); stepclass(Y ~ ., data=DATA, method="lda")
- ロジスティック回帰
-- MASSライブラリ中のstepAIC関数を使ってステップワ...
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])
として列を指定してループ処理する．
* mtext の使いかた 2006 12 28 [#r9850e3f]
「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")
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...
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")
ページ名:
