Site admin: motohiro PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 7.4.33
HTML convert time: 0.069 sec.
個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-m.at.ias.tokushima-u.ac.jp. までご連絡ください．
永田靖著『統計的多重比較法の基礎』サイエンティスト社
bartlett.test(p.37$x ~ p.37$name)
nagata.aov.R 永田靖著『統計的多重比較法の基礎』サイエンティスト社
tapply(p.37$x, list(p.37$name), var)
これは対象群を中心に，それと処理群の平均値の差を調べる
永田 靖 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"))
このサイトを参考にさせていただきました．以下引用です！(石田)．
対数オッズは，通常「ロジット logit」と呼ばれる．
log[ｐ／(1‐ｐ)]＝βｘ＋α
式を書き直すと，
ｐ＝[1＋exp{-(βｘ＋α)}]-1
と表わすことができる．
２群のロジットの差をとると，
logit(ｐ1)‐logit(ｐ2)＝log[ｐ1(1-ｐ2)/[ｐ2(1-ｐ1)]＝β
となる．
これは対数オッズ比log(ψ)がβであることを示している．
すなわち，ψ＝exp(β) である．
Faraway p.27, p.35
Faraway 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)
Faraway 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
---
Faraway 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
Faraway 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)
Faraway p.77
関数 svd の対象は residuals
このサイトを参考にしました．
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.
Faraway p.80
pchisq(deviance(mods), df.residual(mods), lower =F) # [1] 0.003762852 pchisq(deviance(mods), df.residual(mods)) # [1] 0.9962371
< 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
Faraway p.89
data(nes96) # xtabs 使う xtabs(~ PID + educ, nes96) (partyed as.data.frame.table(xtabs( ~ PID + educ, nes96)))
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
Faraway p.98
cutinc さらに ここ も参照
Faraway 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
以下はこのサイトからの引用です！ 石田
帰無仮説どおりのランダムサンプリングによって得られたデータの分布を、 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が同じである仮説同士は同じぐらい有意 であると言うことができる。
Faraway p.30, p.122 - 123
Faraway p.186
iid は独立同型分布(independently and identically distributed)
Faraway p.273
expand.grid(temp = seq(-3, 3, 0.1), ibh = 0, ibt = 0)
ibh, ibt の値は平均(標準化された平均値)に固定し，temp を[-3,3] の間で変化
Wood, GAM, p.5 下 - p.6
例えば平均値が 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))
Wood, p.30
推測されたパラメータが，切片しか含まないモデルから生成された確率。
Wood, GLM, p.33
モデルに不必要な説明変数が多く存在する
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 Analyses Using R p.165
offset の意味 2006 04 28
Wood, GLM, p.45
モデルの含める際，そのパラメータが 1 と仮定する。
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 Using R 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
Everitt & Hothorn, p.68
pairs(means[,-1], panel = function(x, y)
{text(x,y, abbreviate(levels(skulls$epoch)))})
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 p. 831
Everitt & Hothorn, p.98
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個。
Faraway p.41
modl<- glm(cbind(dead, alive) ~ conc,
family = binomial, bliss)
summary(modl)
modl$dev
modl$null
とする．
工学のためのデータサイエンス, 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)
# 高頻度では、データ区間が飛び飛びになっているので、これを補正する
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;
}
}
##
青木先生の説明を参照
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.）。
x <- rnbinom(30, size = 5, prob = 0.5)
ks.test(x,"pnbinom", size = 5, prob = 0.5)
# 間瀬先生著 p.51 下
n<- 4+9+16+13+9+7+5+4+3
x<- rnbinom(n, size = 10.59, prob = 0.7609)
ks.test(x, "pnbinom", size = 10.59, prob = 0.7609)
table(x)
ks.test(table(x), "pnbinom", size = 10.59, prob = 0.7609)
ks.test(observed, "pnbinom", size = my.k,
prob = (my.k / (my.k + mean(x))))
ここで size は ２項分布のパラメータ k (= size) prob は prob<- size/(size + mu) で計算される．
ks.test(rnorm(100), "pnorm") ks.test(rpois(100, 5), "ppois", lambda = 5)
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 p.214, p.511, p.517, p.539 も
Faraway p.121 - 122.
このような記事もあった また次のような記事もあった
From: Dan Kehler (kehler@mathstat.dal.ca) Date: Sat 29 May 2004 - 02:54:48 EST * Next message: Richard Valliant: "[R] dotchart questions" * Previous message: Ted Harding: "RE: [R] Generate a sequence of random integer values" Message-id:<Pine.GSO.3.96.1040528133631.18220B-100000@chase> Using R 1.8.1, and the negative binomial glm implemented in MASS, the default when using anova and a chi-square test is to divide the deviance by the estimated dispersion. Using my UNIX version of S-plus (v 3.4), and the same MASS functions, the deviances are *not* divided by the estimated dispersion. Firstly, I'm wondering if anyone can enlighten about the correct procedure (I thought the F-test was more appropriate when dispersion is estimated)? Secondly, after a bit of muddling with the negative binomial pdf, I concluded that, like for the Poisson, phi is actually 1. This result is borne out by simulations. Is this correct? # an example in R 1.81 with library(MASS) y<-rnegbin(n=100,mu=1,theta=1) x<-1:length(y) # 石田 追加，これは便宜上の説明変数を作ったに過ぎない
model<-glm(y~x,family=neg.bin(1))
summary(model)$dispersion
[1] 1.288926
anova(model,test='Chisq")
#...
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 99 102.038
x 1 0.185 98 101.853 0.705
# But the "real" chi-square probability is
1-pchisq(0.185,1)
[1] 0.6671111
1-pchisq(0.185/1.288926, 1) # 石田 追加
[1] 0.7047963 # 石田 追加
Thanks in advance,
Dan
これは，フルモデルにたいして，ヌルモデルを比べ，観測値の変動をヌルモデルだけで どれだけ説明できるかを示したものである． フルモデルは，本来，deviance がゼロのはずであるので，ヌルモデルとフルモデルの デヴィアンスの差とは，ヌルモデルのデヴィアンスのことではないのか．
また
x<-1:length(y)
は単変量であって，これを使って回帰を行ってもフルモデルを設定したことには ならないのではないか．そもそもフルモデルは自由度 ０ のはずである． その意味で，ヌルモデルそのものがある分布に適しているかどうかを調べるには 先の方法でも十分かも知れない． 下の場合，もともと負の二項分布にのみ従った乱数で，他に変動要因が無い． 尤も，この場合，誤差は負の二項分布の誤差に従っているということになるのか?
y<- 1:length(x)
x.glm.nb2<- glm.nb(x ~ y)
anova(x.glm.nb2, test = "Chi")
これは，以下に等しい
pchisq( x.glm.nb$null.deviance -x.glm.nb$deviance ,
x.glm.nb$df.null - x.glm.nb$df.residual,
lower = F)
anova(x.glm.nb, x.glm.nb2, test = "Chi")
以下の数値とは異なる数値が出る
x.glm.nb2$null.deviance # NULL モデルの deviance
x.glm.nb$df.null # NULL モデルの自由度
pchisq( x.glm.nb2$null.deviance , x.glm.nb$df.null,
lower = F)
Faraway p. 124
response, pearson, working とあるが， 解析結果の $residuals で出力されるのは working その他の residual を得るには関数 residual を使う．
suppressWarnings 関数の引数として，命令を実行する．
Faraway p.29 -30. またカイ自乗分布は，平均 ＝ 自由度，分散 ＝ 2*自由度に従う．
Faraway p.33
Faraway p.36
描画を行う
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))
正確に行うには 青木先生のサイト
http://aoki2.si.gunma-u.ac.jp/lecture/GoodnessOfFitness/normaldist.html 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.ob) )
sum(log.V1.exp)
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
################ リストの中身を全て取り出す
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(bun.p >= 4){
ex.freq<- 5
}else{
ex.freq<- 1
}
は良いが
if(bun.p >= 4){
ex.freq<- 5
}else{
ex.freq<- 1
}
はエラーとなる．
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,], 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)
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
## 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)))
渡部『ベイズ統計学入門』 p.72
モデル分布がベルヌーイ分布である母集団から n = 3 の標本を抽出したとき，その観測値は
c(1,0,1)
であった．母数の最尤推定値を求めよ．
f <- function(x) x^2 * (1-x) optimize(f,c(0.001, 0.999),maximum=T)
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
# 密度関数を独自定義
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: # 関連警告
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
matsuda.121<- array(c(32,86, 11, 35, 61,73,41,70),
dim = c(2,2,2),
dimnames = list(height=c("over","lower"),
diam = c("lower","heigher"),
species = c("sagrei", "distichus")))
class( matsuda.121 )<- "table"
matsuda.121
# p.129
matsuda.121.loglm <-
loglm(~ height * diam * species, data = matsuda.121 )
stepAIC(matsuda.121.loglm )
Step: AIC= 14.03
~height + diam + species + height:species + diam:species
Df AIC
<none> 14.026
- height:species 1 22.431
- diam:species 1 24.632
Call:
loglm(formula = ~height + diam + species +
height:species + diam:species,
data = matsuda.121, evaluate = FALSE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 2.025647 2 0.3631921
Pearson 2.017364 2 0.3646994
同じことだが
matsuda.121.loglin <- loglin( matsuda.121, list(c(1,3), c(2,3)) ) matsuda.121.loglin
c(1,2) を外してよいか．つまり heigh と diam は独立だとみなし，無視してよいか 実際に外した結果を，飽和モデルと比較する 0.3631921 なので，外してよい ↓
1 - pchisq(matsuda.121.loglin$lrt, matsuda.121.loglin$df) 1 - pchisq(matsuda.121.loglin$pearson, matsuda.121.loglin$df)#
よって species さえ分かれば，height や diam の間に特に関係は無い．
次のように，デフォルト以外のライブラリの場所を確認
Sys.getenv("R_LIBS")
必要があれば
Sys.setenv(R_LIBS="C:/R")
install.packages("MyPackages", lib = "C:/R")
# 第二引数はオプション
library("MyPackages")
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
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("日本語ひらがら"))
例えばEveritt & Hothorn, A Handbook of Statistical Analyses Using R 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)])
# 相関行列 # 条件指定を変更
Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.75 下
Everitt & Hothorn, A Handbook of Statistical Analyses Using R 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)]
Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.51
cmh_test( classification ~ treatment, data = Lanza,
scores = list(classification = c(0, 1, 6, 17, 30)),
subset = Lanza$study == "I")
cmh_test( classification ~ treatment, data = Lanza,
scores = list(classification = c(0, 1, 6, 17, 30)),
subset = Lanza$study == "II")
# p.82
psymb<- as.numeric(clouds$seeding)
plot(rainfall ~ cloudcover , data = clouds, pch = psymb)
abline(lm(rainfall ~ cloudcover, data = clouds,
subset = seeding == "no"))
abline(lm(rainfall ~ cloudcover, data = clouds,
subset = seeding == "yes"), lty = 2)
legend("topright" , legend = c("No seeding", "Seeding"),
pch = 1:2, lty = 1:2, bty = "n")
Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.106
また overdispersion が生じるのは，しばしばデータが独立でないため
Everitt & Hothorn, A Handbook of Statistical Analyses Using R p.117
Everitt & Hothorn, A Handbook of Statistical Analyses Using R p. 166
例えば以下のコードは
BtheB[, grep("bdi", names(BtheB))]
BtheB データフレームで，列名に "bdi" がついている列番号を抽出している．
さらに以下は
subset(BtheB, treatment == "TAU")
[, grep("bdi", names(BtheB))]
BtheB データで，列 treatment が "TAU" であるサブセットから， "bdi" を名前として含む列を取り出している．
Everitt & Hothorn, A Handbook of Statistical Analyses Using R chapter 10 p. 165
data("BtheB", package = "HSAUR")
BtheB$subject<- factor(rownames(BtheB))
nobs<- nrow(BtheB)
BtheB.long<- reshape(BtheB, idvar = "subject",
varying = c("bdi.2m", "bdi.4m", "bdi.6m", "bdi.8m"),
direction = "long")
# 石村貞夫「分散分析のはなし」の例．データは p.77
# データフレームの作成
# kaeru<- read.csv("africanFlog.csv",header=F)
# 始めに横長形式のデータとして与えられていたとする．
frogs<- data.frame(stage = c("A1", "A2", "A3", "A4", "A5"),
sample1 = c(12.2, 22.2, 20.8, 26.4, 24.5),
sample2 = c(18.8, 20.5, 19.5, 32.6, 21.2),
sample3 = c(18.2, 14.6, 26.3, 31.3, 22.4))
frog.res<- reshape(frogs, idvar = "stage",
varying = list(names(frogs[2:4])),v.names = "data",
direction = "long")
# ここで作成される time 変数は，各数値の順番を表す
# (時間ごとにデータを取っている場合などに重要)
frog2<- frog.res[order(frog.res$stage),] # 並び替え frog2
以下は Rjp Wiki よりの引用です！ 石田
判別分析では、モデルに含まれるパラメーター数が多くなればなるほど、あてはめ誤差（残差平方和）は小さくなるということがあります。そこで、単に残差平方和の大小を比較するだけでなく、パラメータ数も考慮したモデルの選択を行います。その際の基準となるものがAICです。wleライブラリ中のmle.aic関数を使ってモデルを選択します。
library(MASS) library(wle) result<-mle.aic(lda(class~.,data=DATA)) summary(result,num=20)
mle.cvでやった場合にも判別率は出力されませんでした。AIC値の代わりにcv値が出力されていますが、どういう値なのかよくわかりません。 上記に続いて、判別関数を求める手順、交差妥当性（クロス・バリデーション）のチェックをやって判別率を比べる手順を示します。
(z<-lda(class~.,data=DATA)) apply(z$means%*%Z$scaling,2,mean)
この結果、それぞれの変数の係数と定数が出力されます。
predict(z)$x predict(z)$class predict(z)$posterior
と入力すればこの判別関数を用いた判別得点、各個体が判別されたグループのラベル、各個体がどのグループに判別されているかに関する事後確率（0から1）が出力されます。判別率を計算するためのクロス表の出力は次のようにします。
table(data[,],predict(z)$class)
交差妥当性のチェックを行ったクロス表は
DATA.CV<-lda(class~.,data=DATA,CV=T) (lda.tab<-table(DATA[,],data=DATA.CV$class))
で出力できます。判別率、誤判別率はそれぞれ
sum(lda.tab[row(lda.tab)==col(lda.tab)])/sum(lda.tab) sum(lda.tab[row(lda.tab)!=col(lda.tab)])/sum(lda.tab)
交差確認（cross validation）はデータを2分し,一方のサンプルで判別式をつくりその判別式で残りのサンプルを判別して, 判別式の有効性を見る方法です。データを2分する以外に3,4,5,10分する方法もありますが、どれを用いるかは明確な基準がなくデータサイズに依存するそうです。 この特殊なケースとして1 つ取っておき法（reave-one-out cross-validation）があります。 データから1つの個体を取り除いて判別分析を行い、取り除かれた個体で判別モデルの評価を行うという作業を全ての個体に繰り返して行います。 「R」でこれが出来るということなので紹介します。ここでは「R」の中にあるirisデータを例に使います。
library(MASS) iris.CV<-lda(Species~.,data=iris,CV=T) (lda.tab<-table(iris[,5],iris.CV$class))
デフォルトではCV=FALSEとなっているのでCV=Tにすると1つ取っておき法（reave-one- out cross-validation）ができます。 ちなみに判別関数lda(Linear Discriminant Analysis),qda(Quadratic Discriminant Analysis)には1つ取っておき法（reave-one-out cross-validation）以外の交差確認機能は用意されていないそうです。 詳しい説明は,同志社大学Jin先生のページhttp://www1.doshisha.ac.jp/~mjin/R/の「Rと判別分析」を見ると書いてあります。
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)
mypredict.ldaはclass labelだけを返させるのに必要な手続きらしい。
library(klaR); stepclass(Y ~ ., data=DATA, method="lda")
library(MASS); stepAIC(glm(Y ~ ., data=DATA,
family=binomial))
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<- 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])
として列を指定してループ処理する．
「Rの基礎とプログラミング技法」（シュプリンガージャパン）169ページのグラフィックス
ligges.R
par(oma=rep(3, 4), bg="white")
plot(c(0, 1), c(0, 1), type="n", ann=FALSE, axes=FALSE)
rect(-1, -1, 2, 2, col="white")
box("figure")
par(xpd=FALSE)
rect(-1, -1, 2, 2, col="white")
box("plot", lty="dashed")
text(.5, .5, "plot region", cex = 1.6)
mtext("figure region", side=3, line=2, adj = 1, cex = 1.4)
for (i in 1:4)
mtext(paste("inner margin", i), side=i,
line=1, outer= FALSE)
for (i in 1:4)
mtext(paste("outer margin", i), side=i,
line=1, outer=TRUE)
mtext("device region", side=3, line=2,
outer=TRUE, adj = 1, cex = 1.2)
box("outer", col="black")