Site admin: motohiro PukiWiki 1.5.1 © 2001-2016 PukiWiki Development Team. Powered by PHP 7.4.33
HTML convert time: 0.012 sec.
#content
* [[ ]] の使い方 [#z6109ecf]
例えば以下のようなデータで ([[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.200)
> z
0.8 1 1.2 1.4 1.6 2.5
0 6 4 2 2 0 0
1 1 1 4 4 4 2
>
> z[[1]] // 一列目の一行目
[1] 6
> z[[2]] // 一列目の二行目
[1] 1
> z[[3]] // 二列目の一行目
[1] 4
> z[[4]] // 二列目の二行目
[1] 1
> z[[5]] // 三列目の一行目
[1] 2
> z[[6]] // 三列目の二行目
[1] 4
>
> z
このとき以下は,行名を
> dimnames(z)[[1]]
[1] "0" "1"
次は列名を表示する.
> dimnames(z)[[2]]
[1] "0.8" "1" "1.2" "1.4" "1.6" "2.5"
今は二次元のデータだから,三次元を指定するとエラーとなる.
> dimnames(z)[[3]]
Error in dimnames(z)[[3]] : subscript out of bounds
* boxplot [#qb7738a5]
例えばカッコーの巣であれば，
> source("/source/statistics/R/book/DAAG/data/cuckoos.R")
> cuckoos
length breadth species id
1 21.7 16.1 meadow.pipit 21
...
120 21.0 16.0 wren 238
> boxplot(length ~ species, data=cuckoos)
>
ちなみに，外れ値はレンジで指定するが，デフォルトでは range = 1.5 となっている。詳細は『データサイエンス入門』p.14
またrug() で，y軸右(side=4)ないし左(side=2)に，プロットに対応したバーコードを表示できる。
ラグプロットについてはここを参照
http://www.okada.jp.org/RWiki/index.php?%A5%E9%A5%B0%A5%D7%A5%ED%A5%C3%A5%C8
なお boxplot と split, rug を併用した例は [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.225
* リストの先頭要素を取る [#r796e7ac]
> x > x[1:2]
$string
[1] "a"
$integer
[1] 1 2
> x["string"]
$string
[1] "a"
> x[-1]
$integer
[1] 1 2
$float
[1] 3.14
* qt下側確率 [#uf69c824]
構文は qt(下側確率，自由度) 坪田 p.29
例えば qt(0.05,10) とすると，qt()は分位点関数なので，正規曲線の総面積の累積が5%となる x の座標点
> qt(0.05,10)
[1] -1.812461
を与える。上側は
> qt(0.95,10)
[1] 1.812461
* Excel データの読み込み [#va762f1e]
RODBCパッケージを利用する。
以下は，excelファイルを読み込むだけでなく，更新なども行なう。
library(RODBC)
connExcel
このあと
data
などとクエリーを行なえる。使い終われば
odbcClose(connExcel)
* 新しいグラフィックデバイス [#t51cebd3]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] , p.22,35
X11() を実行。サイズを指定するには X11(width=7.5,height=4)
現在アクティヴなグラフィックデバイスを閉じるには
dev.off()
特定のデバイス, 例えば ３ 番を閉じるには
dev.off(3)
デバイスの変更, 例えば ２ 番に切替えるなら
dev.set(2)
* 特定の行の取り出し [#l372d554]
例えば [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] , p.179 の砂糖の平均を，コントロール群だけ取り出すなら
> sugar
weight trt
1 82.0 Control
2 97.8 Control
3 69.9 Control
4 58.3 A
...
12 48.9 C
> mean(sugar$weight[sugar$trt == "Control"])
[1] 83.23333
>
あるいは
> splitplot[1,]
Block Insect Mollusc Rabbit Lime Competition Nutrient Biomass
1 A Sprayed Slugs Fenced Unlimed Control N 7.95395
...
のようなデータフレームがあるとして，Block が A の項目は
splitplot[splitplot$Block=="A",]
として得られる。
~/mysplus/crawley.s , p.352
次の例は，term.cos で，行名が Berg である行の各列から，0.999 以上の項目があれば取り出す
berg >= 0.999]
# あるマトリックスの，ある名前の行の (1列目の) 数値を見る
stif.keller.without.svd.Doc3[rownames(stif.keller.without.svd.Doc3)=="Gouverneur", 1]
* 図における曲線，座標設定などの詳細 [#j0fca1af]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.182
なお plot で引数に type="n"を指定すると，点は刻印されない。これに続いて
points() を実行する。
二つのカテゴリについて，変数 Y と 変数 X の対応をプロットする。[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.134
* aggregate の使い方 [#cbfa9ae1]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.43 例えば
yield block shade plot
101 north none north.none
108 north none noth.none
のようなデータがあった時，ブロックとシェードをまとめた平均を作成する場合
aggregate(yield, by=list(block,shade), mean)
block shade meanOFyield
east none 99
他に Maindonals p.225
* rank 落ち [#s5b1c832]
// 例えば一つの水準に三つの項目があり，どれか一つだけを選択する場合，
// 任意に選んだ二つの水準のそれぞれについて，どちらかが選択されているかどうかが分かれば，
// 残り一つの水準について選択されているかどうかは自動的に分かる。
//
// つまり一つの水準の値は，残り２水準の値によって計算される。
// すなわちある水準のベクトルは，残り２水準のベクトルの和で表され，
// それらに従属している。
* ニューラルネットワークその２ [#l4142b36]
豊田「非線形多変量解析」87ページを実行するには
data(iris)
library(nnet)
iris.2 # このままだと元のデータフレムの情報として，Species 変数のレベルが三つになっている。
levels(iris.2$Species)
#[1] "setosa" "versicolor" "virginica"
そこで三つ目のレベルを欠損値とする。これによりレベルは二つとして扱われる。
levels(iris.2$Species)[3]
iris.nnet iris.nnet.pre table(iris.2$Species, iris.nnet.pre)
summary(iris.nnet)
* 決定木でラベルを表示する [#hdf800ef]
text 関数が有効な時と，有効でない時がある。これは，データフレームの作成方法に問題がある。
有効
senkei senkei.tree plot(senkei.tree,type="u")
text(senkei.tree)
summary(senkei.tree)
無効
base base1 <- base.tree + man + mit + nicht
+noch+nun+nur+oder+schon+sich+so+und+zu, data=base1)
plot(base.tree,type="u")
text(base.tree)
summary(base.tree)
ところが，一度
text(base.tree,labels=col.names(base.tree))
Error in text(xy.coords(x, y, recycle = TRUE), labels, adj, pos, offset, :
invalid pos value
In addition: Warning message:
NAs introduced by coercion
を実行 (エラーで無表示となる)後，再度
text(base.tree)
を行なうと，ラベルが付与された。
* 多重共線性の発見 [#p762e5ab]
朝野 p.108より
説明変数の相関行列 R について R の固有値を求める。このとき固有値の積は，R の行列式に一致する。偏回帰係数の分散は，固有値の逆数の和と一致する。
* 自己組織化マップ [#pffef85e]
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=R%A4%C7SOM%28%BC%AB%B8%CA%C1%C8%BF%A5%B2%BD%A5%DE%A5%C3%A5%D7%29&word=som
#
library(som)
library(class)
ex gr#ex[1,]
testplot(test)
# 六角形配置にしたがって円を配置する (下記の 12*8は xdim, ydimに対応)
symbols(test$grid$pts[,1],test$grid$pts[,2],circles=rep(0.5,8*12),inches=F,add=T)
# knn1を使って ex 中のデータが SOMで割り当てた test$code (8*12ある)のどれに近いのかを調査
bins # exから抜粋した値がどこに割り当てられるのかを示す。
# 全てをプロットすると収集がつかなくなるので seq(101, 1000, by=25)だけをプロットしてみる
text(test$grid$pts[bins[seq(101,1000,by=25)],]
+ rnorm(36,0,0.2),col="blue",as.character(seq(101,1000,by=25) ))
* データフレームの操作 [#mf1618ef]
stifstom4[stifstom4
stifstom4[,stifstom4[,c(4:16)]
stifstom4[,c(4:16)][stifstom4[,c(4:16)] #stifstom4[,c(4:16)][stifstom4
stifstom4[stifstom4 > 1]
* Dalgaard の本のデータ解析 [#c1ad686c]
research/statistics/Intro.R.R
* matrixの作り方 [#n6b7dd23]
例えば自然科学の統計学 p.38 の最小二乗法の計算
x A B X t(X) %*% X
* 行列式を使った方程式の解 [#l7fd71d4]
自然科学の統計学 p.67
3x - 2y = 7
x + 3y = -5
は R で以下のように解く
x1 X #X = matrix(c(3,-2,1,3),ncol=2,byrow=T) でももちろん良い。
Y = c(7,-5);Y
solve(X) %*% Y
# あるいは
solve(X,Y)
* 松原望 蟻が黒い確率 [#g145d398]
松原望『実践としての統計学』 p.16, p.163, p.222
B1 = 100% 蟻は黒い, とする。
B2 = 99 % の蟻は黒い, とする。そしてそれぞれの確率を
p(B1) = p(B2) = 0.5
すると
p(D|B1) = 1, 蟻は100%黒いという仮定が成立して，その時に黒い蟻を見つける確率 。
p(D|B2) = 0.99, 蟻は99%黒いという仮定が成立して，その時に黒い蟻を見つける確率。
このとき p(B1|D) , 黒い蟻を見つけた時に，B1 が正しい確率は，
p(B1|D) = p(B1 * D) / p(D)
= p(B1) * p(D|B1) / p(B1) * p(D|B1) + p(B2) * p(D|B2)
= 0.5 * 1 / ( 0.5 * 1 + 0.5 * 0.99)
と計算される。
* 作業ディレクトリにあるすべてのオブジェクトの削除 [#u8b09aa4]
remove(objects()) # S-Plus
rm(list = ls()) # R
* 図のサイズ調整 [#ve59356c]
# dev.off()
# plot.new() 場合によっては描画前に必要となる
par.old par.old plot(x,y)
par(par.old)
Dalgaard p. 28, 30, 71
Venables p. 84
これに外周を加えるのが
oma
* ブラウザからヘルプを見る [#je7d86ba]
help.start()
* R のライブラリ類の場所 [#f343fc64]
/usr/local/lib/R/library/
* データテーブルを並び替える stab [#u850d32b]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.37 以下のようなデータがあった場合
> jobs
BC Alberta Prairies Ontario Quebec Atlantic Date
1 1752 1366 982 5239 3196 947 95.00000
2 1737 1369 981 5233 3205 946 95.08333
3 1765 1380 984 5212 3191 954 95.16667
...
Jobs
によって，次のようなデータに変更できる。
> Jobs
values ind
1 1752 BC
2 1737 BC
3 1765 BC
...
142 944 Atlantic
143 943 Atlantic
144 942 Atlantic
* abline と spline, lines の違い [#p65d45d0]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] Chapter3, p.53, p.54, p.55
abline(x,y) は切片が x 傾きが y の直線を引く
splineは lines と組み合わせて，曲線
lines(spline(x,y))
他に
lines(lowess(x,y),) p.35
も使われる。
* pretty, axis の使い方 [#n74c2833]
pretty は，第一引数をもとに，第二引数で指定された数のベクトルを作成する。
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.59, p.117
* set.seed [#peb6a616]
set.seed(x) は，R起動には，x をランダムな値に設定して生成される乱数のリストである。
したがって，同じ x の値で実行すれば，同じ乱数のセットが得られるので，
結果として，同じ結果が得られる。[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.64
* split [#s1d4876e]
第一引数で指定したベクトルを，第二引数で指定したカテゴリーで，リストに分ける。
Maindpnald p. 134, p.213
* text関数の使い方 [#a0c82837]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.134 より
plot(weight.split$hb ~ volume.split$hb,
pch=16, xlim=range(volume), ylim = range(weight),
ylab="Weight (g)", xlab="Volume (cc)")
text(weight.split$hb ~ volume.split$hb, labels=c(1:7), adj=c(1.5,-1.5))
ここで adj は(x,y)点に対して x がマイナスなら，左に，また y がマイナスなら上にテキストが付記される。
* 回帰分析で，変数間の相関係数を出すオプション [#f7ea1c56]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] # p.137
allbacks0.lm summary(allbacks.lm, corr=T)
また個々のサンプルを抜いた場合の影響度を調べるには [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.139
round(coef(allbacks.lm0), 2)
round(lm.influence(allbacks.lm0)$coef, 2)
* 変数間の線型性を求める alias [#j5f121ed]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.166
* 分散分析，回帰分析で base ラインを変更する relevel [#sc3be186]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.90 p.178 p.210
* option の使い方 [#e595e701]
Maindonals p.178
* 散布図の外にプロットを打つ points(, ..., xpd =T) [#cf96f320]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.202
* nlme マルチレベルでの固定効果と，ランダム効果 [#w367daae]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.226
science.lme summary(science.lme)$tTable
* カイ二乗検定における有意確率(両側) [#teb7b8ff]
青木先生のサイト
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/one-two.html
* 浜田知久馬 ロジスティック回帰分析 [#rebda9b4]
hamada.R を参照 p.141
x y z # z
# x y
# [1,] 95 5
# [2,] 90 10
#
fac # fac.z fac.z gl(2,1,2,fac); fac.z
# fac.z
# [1] A B
# Levels: A B
test.glm glm(z ~ fac.z, family=binomial)
summary(test.glm)
summary(test.glm)$coef[1,3]^2 # χ二乗
summary(test.glm)$coef[2,3]^2 # χ二乗
#exp(abs(summary(test.glm)$coef[2,1]))
exp(summary(test.glm)$coef[2,1])
anova(test.glm, test="Chisq")
fisher.test(z)
# p.143
x y z #> z
# x y
#[1,] 5 95
#[2,] 2 48
#[3,] 8 42
fac fac.z #> fac.z
#[1] A B C
#Levels: A B C
test.glm summary(test.glm)
summary(test.glm)$coef[1,3]^2 # χ二乗
summary(test.glm)$coef[2,3]^2 # χ二乗
summary(test.glm)$coef[3,3]^2 # χ二乗
exp(summary(test.glm)$coef[2,1])
exp(summary(test.glm)$coef[3,1])
anova(test.glm, test="Chisq")
# p.144
x y z
fac
test.glm summary(test.glm)
summary(test.glm)$coef[1,3]^2 # χ二乗
summary(test.glm)$coef[2,3]^2 # χ二乗
anova(test.glm, test="Chisq")
* array の作り方 [#p6acf9ec]
Meyer English Corpus Linguistics p.134
Type Form 1-4 Words 5 or more Words
PT Simple NP 216 23
PT Gen. NP 0 0
PT Post.Mod. 1 7
APOS Simple NP 52 10
APOS Gen. NP 18 9
APOS Post.Mod. 12 97
Meyer array(c(
216, 0, 1,
23, 0, 7,
52, 18, 12,
10, 9, 97),
dim = c(3,2, 2),
dimnames =
list(
Form = c("SimpleNP", "GenNP", "PostMed"), # = c(1)# = c(1)
Length = c("lessWords","moreWords"), # = c(2)
Type= c("PT", "AP") # = c(3)
))
class(Meyer) ftable(Meyer)
c(3, 2, 2) で 3 行 2 列 2 層，入力順に，層(行→列)→層(行→列)と埋める
> Meyer
, , Type = PT
Length
Form lessWords moreWords
SimpleNP 216 23
GenNP 0 0
PostMed 1 7
, , Type = AP
Length
Form lessWords moreWords
SimpleNP 52 10
GenNP 18 9
PostMed 12 97
> ftable(Meyer)
Type PT AP
Form Length
SimpleNP lessWords 216 52
moreWords 23 10
GenNP lessWords 0 18
moreWords 0 9
PostMed lessWords 1 12
moreWords 7 97
>
####################
Meyer array(c(
216, 23,
0,0,
1,7,
52,10,
18,9,
12,97),
dim = c(2,3, 2),
dimnames =
list(
Length = c("Less", "More"),
Form = c("SNP", "GNP","PMD"),
Style= c("PT", "AP")
))
c(2, 3, 2) で２行３列，入力順に行→列と埋める
2*3 を埋め尽くすと，次に次元に移り，2 * 3 を再度埋めていく
> Meyer
, , Style = PT
Form
Length SNP GNP PMD
Less 216 0 1
More 23 0 7
, , Style = AP
Form
Length SNP GNP PMD
Less 52 18 12
More 10 9 97
class(Meyer) # "array"
class(Meyer)<- "table"
class(Meyer) # "table"
アレイをテーブルに変更する
> ftable(Meyer)
Style PT AP
Length Form
Less SNP 216 52
GNP 0 18
PMD 1 12
More SNP 23 10
GNP 0 9
PMD 7 97
* ログリニア分析 [#q5ce84c6]
始めに[[ここ>#matsuda121]]を参照．[[ここ>#loglin2]と[[ここ>#loglm]]も．[[ここ>#social.mobility2]も．さらにMeyer.Rを参照。
Meyer.Howa
ここで list(c(1,2,3)) は list(1:2:3)としても同じで，三次交互作用までを含んだモデル。
list(c(1,3), 2) ならば，1,3 相互の交互作用を含めるが，
2 については，1,3 とは独立の因子とみなす。
list(1,2,3) とすれば，すべてが独立のモデルである。
なお階層モデルで三次の交互作用までを含むということは，2 次，また単独の効果を含む。柳井/緒方『統計学-基礎と応用』,p.150
飽和モデルの結果をみる
# meyer p.135 表，３次の効果を含む飽和モデル，自由度０で説明力は無い
Meyer.Howa$lrt # 尤度比
1-pchisq(Meyer.Howa$lrt,Meyer.Howa$df) # 有意度
Meyer.Howa$pearson# χ二乗値
1-pchisq(Meyer.Howa$pearson,Meyer.Howa$df) # 有意度
Meyer.Howa$df # 自由度
Meyer.Howa$margin # 選択されたモデル。この場合は飽和モデル
#[[1]]
#[1] "Form" "Length" "Type"
###################
# meyer p.135 表，３次の効果を除いた場合
# K = 2
Meyer.123<- loglin(Meyer, list(c(1,2),c(2,3),c(1,3)), param=T)
#c(1,2,3) を除いたモデルの適合度は 0.926 ，つまり有意でないので，３次の交互作用を除いたモデルは適合している。
すなわち３次の適合度は省いて構わない。
Meyer.123$lrt #0.1530465
1-pchisq(Meyer.123$lrt,Meyer.123$df) #0.9263314
1-pchisq(Meyer.123$pearson,Meyer.123$df)
Meyer.123$pearson
Meyer.123$df # 2
Meyer.123$margin
#[[1]]
#[1] "Form" "Length"
#[[2]]
#[1] "Length" "Type"
#[[3]]
#[1] "Form" "Type"
############ meyer p.135 頁の表, ３次に加えて，さらに２次を全て省いた場合
# K = 2
Meyer.1.2.3<- loglin(Meyer, list(1,2,3), param=T)
Meyer.1.2.3$lrt # 488.0099
1-pchisq(Meyer.1.2.3$lrt,Meyer.1.2.3$df) # 0
1-pchisq(Meyer.1.2.3$pearson,Meyer.1.2.3$df)
Meyer.1.2.3$pearson
Meyer.1.2.3$df # 7
Meyer.1.2.3$margin
#[[1]]
#[1] "Form"
#[[2]]
#[1] "Length"
#[[3]]
#[1] "Type"
* loglm によるログリニア，対数線形分析 [#ge2b6b48]
&aname(loglm);
始めに[[ここ>#matsuda121]]を参照．[[ここ>#loglin2]]も．[[ここ>#social.mobility2]]も
loglm 使うためにMASSパッケージを取り込み
library(MASS)
初めに飽和モデルを分析する。
Meyer.loglm<- loglm(~Length*Form*Type, data=Meyer)
次に，AIC が最小になるモデルを選び出す。
# 飽和モデルから，効果を減らしてく
Meyer.loglm<- loglm(~Length*Form*Type, data=Meyer)
anova(Meyer.loglm)
Meyer.loglm2<- update(Meyer.loglm, .~. -Form:Type)
anova(Meyer.loglm2)
stepAIC(Meyer.loglm, ~.^2, data=Meyer)
#######
Start: AIC= 24
~Length * Form * Type
Df AIC
- Length:Form:Type 2 20.153# Length:Form:Type を省いた方が小さい
<none> 24.000
Step: AIC= 20.15
~Length + Form + Type + Length:Form + Length:Type + Form:Type
Df AIC
- Length:Type 1 19.978# Length:Type を省いた方が小さい
<none> 20.153
+ Length:Form:Type 2 24.000
- Length:Form 2 145.147
- Form:Type 2 153.044
Step: AIC= 19.98
~Length + Form + Type + Length:Form + Form:Type
Df AIC
<none> 19.978 # これ以上 AIC は小さくならない
+ Length:Type 1 20.153
- Length:Form 2 255.045
- Form:Type 2 262.943
Call:
loglm(formula = ~Length + Form + Type + Length:Form + Form:Type,
data = Meyer, evaluate = FALSE)
Statistics:
X^2 df P(> X^2)
Likelihood Ratio 1.977911 3 0.577004
Pearson NaN 3 NaN
>
#####x よって最適なモデルは
# ~Length + Form + Type + Length:Form + Form:Type
#####
別データによる分析
それぞれの効果の推定値には，ダミーコーディングによる場合と，ANOVAコーディングによる場合がある。
rocket<- array(data=c(5,7,8,9,3,21,7,9,6), dim=c(3,3))
dimnames(rocket)<- list(c("A1","A2","A3"),c("B1","B2","B3"))
rocket
# モデル選択
model0<- loglin(rocket, list(c(1, 2)), param=TRUE)
model1<- loglin(rocket, list(1, 2), param=TRUE)
model0
model1
p0<- 1-pchisq(model0$lrt, model0$df)
p1<- 1-pchisq(model1$lrt, model1$df)
p0
p1
AIC0<- model0$pearson-2*model0$df
AIC1<- model1$pearson-2*model1$df
AIC0
AIC1
# ここでは飽和モデルを採用し分析してみる
rocket.df<- as.data.frame.table(rocket)
colnames(rocket.df)<- c("kozoku","yoko","dosu")
> loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
2 iterations: deviation 0
$lrt
[1] 0
$pearson
[1] 0
$df
[1] 0
$margin
$margin[[1]]
[1] 1 2
$margin[[2]]
[1] 1
$margin[[3]]
[1] 2
$param
$param$"(Intercept)"
[1] 1.990005
$param$"1"
A1 A2 A3
-0.07248058 -0.24275578 0.31523636
$param$"2"
B1 B2 B3
-0.11174159 0.12344831 -0.01170672
$param$"1.2"
B1 B2 B3
A1 -0.1963447 0.1562521 0.04009266
A2 0.3104027 -0.7720850 0.46168230
A3 -0.1140580 0.6158330 -0.50177496
# ダミーコーディングによる効果の推定値
glm0<- glm(dosu ~ kozoku * yoko, data = rocket.df, family = poisson)
anova(glm0, test = "F")
summary(glm0)
#########
> summary(glm0)
Call:
glm(formula = dosu ~ kozoku * yoko, family = poisson, data = rocket.df)
Deviance Residuals:
[1] 0 0 0 0 0 0 0 0 0
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.60944 0.44721 3.599 0.000320 ***
kozokuA2 0.33647 0.58554 0.575 0.565538
kozokuA3 0.47000 0.57009 0.824 0.409689
yokoB2 0.58779 0.55777 1.054 0.291970
yokoB3 0.33647 0.58554 0.575 0.565538
kozokuA2:yokoB2 -1.43508 0.88730 -1.617 0.105800
kozokuA3:yokoB2 0.37729 0.69551 0.542 0.587492
kozokuA2:yokoB3 -0.08516 0.77254 -0.110 0.912227
kozokuA3:yokoB3 -0.62415 0.79657 -0.784 0.433303
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 2.0874e+01 on 8 degrees of freedom
Residual deviance: -4.4409e-16 on 0 degrees of freedom
AIC: 52.681
Number of Fisher Scoring iterations: 3
>
#######
# ANOVA コーディング効果の推定値
loglm(dosu ~ kozoku * yoko, data = rocket.df, family = poisson)$param
loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
#######
> loglm(dosu ~ kozoku * yoko, data = rocket.df, family = poisson)$param
$"(Intercept)"
[1] 1.990005
$kozoku
A1 A2 A3
-0.07248058 -0.24275578 0.31523636
$yoko
B1 B2 B3
-0.11174159 0.12344831 -0.01170672
$kozoku.yoko
yoko
kozoku B1 B2 B3
A1 -0.1963447 0.1562521 0.04009266
A2 0.3104027 -0.7720850 0.46168230
A3 -0.1140580 0.6158330 -0.50177496
########
> loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
2 iterations: deviation 0
$lrt
[1] 0
$pearson
[1] 0
$df
[1] 0
$margin
$margin[[1]]
[1] 1 2
$margin[[2]]
[1] 1
$margin[[3]]
[1] 2
$param
$param$"(Intercept)"
[1] 1.990005
$param$"1"
A1 A2 A3
-0.07248058 -0.24275578 0.31523636
$param$"2"
B1 B2 B3
-0.11174159 0.12344831 -0.01170672
$param$"1.2"
B1 B2 B3
A1 -0.1963447 0.1562521 0.04009266
A2 0.3104027 -0.7720850 0.46168230
A3 -0.1140580 0.6158330 -0.50177496
なおこのデータ，またコーディングについては
http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03-01.pdf
の表3.2「ロケットの発射試験」
この例題の飽和モデルによる推定値は
http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03-04.pdf
の表3.11にダミーコーディング、表3.12にANOVAコーディング
なお SPSS で出力される Partial association は「標準化係数（標準効果）」のこと。
http://www.okada.jp.org/RWiki/index.php?
cmd=read&page=%A3%D1%A1%F5%A3%C1%20%28%BD%E9%B5%E9%BC%D4%A5%B3%A1%BC%A5%B9%29&word=loglin
* LaTeX でテーブル 1 [#h6be20ef]
library(xtable)
x<- matrix(1:24,4,6)
mytable<- xtable(x)
print(mytable)
* LaTeX でテーブル 2 [#s319ea2e]
http://www.okada.jp.org/RWiki/index.php?cmd=read&page=LaTeX&word=TeX%20#content_1_0
* 行列のランクを求める。 [#s806caeb]
qr(X)$rank # 舟尾 p.61
* ベクトルのノルム [#b92c3db0]
渋谷 柴田 p.37
norm if(p=="inf") return (max(abs(x)))
sum(abs(x)^p)^(1/p)
}
これは sqrt(t(x) %*% x) に等しい。
* read.csv 読み込み時に変数名を付加する [#o54b5cb7]
例えば，以下のようにする場合
#x<- read.csv("/research/corpus/file/anborgia.nr.csv",header=F,row.names=1);x
# colnames(x)<- c("line");x
以下のように実行しても構わない。ただし元データが２列で，その一列目を行名と指定している。
その関係で，col.namesの引数は二ついる。
x<- read.csv("/home/ishida/research/corpus/file/anborgia.nr.csv",
header=F,row.names=1,col.names=c("","line1"));x
* 生データから頻度表を作成し，適合度の検定を行なう [#zaf0a874]
まず
> jena > jena
freq
1 15
2 10
3 21
4 20
5 30
...
> kanzel > kanzel
freq
1 19
2 38
3 4
4 23
5 29
...
table( factor(cut(buntyo$V1, breaks=c(0, 5, 10, 15, 25, 35, 40))) )
などとする．
ここでは 0 より大きく 5 以下，5 より大きく10以下，などを区間としている．
としてデータを読み込み，これに頻度カテゴリを表す列を追加する。
tb jena$fac
> jena
freq fac
1 15 (11,16]
2 10 (6,11]
3 21 (16,22]
4 20 (16,22]
5 30 (22,31]
...
> kanzel
freq fac
1 19 (16,22]
2 38 (31,92]
3 4 (0,6]
4 23 (22,31]
5 29 (22,31]
ここで table(kanzel$jena)とすれば，各カテゴリの頻度が求める。
さて，ここでこの二つの頻度について適合度の検定を行なうには，
test chisq.test(test)
とする。
* 図に凡例を付記する [#wa609321]
legend(x軸, y軸, ラベルの指定, fill= 色の指定, pch = シンボル | lty = 線種)
labs cols mypch
legend(7, 10, labs, pch = mypch, col=cols)
x,y 軸は locator(1)として調べれば良い。あるいは引数として与える。
[[John Fox:http://socserv.mcmaster.ca/jfox/]], p.138
legend(locator(1), c('high', 'low'), lty = c(1,3), lwd =c(3,3), pch = c(19,1))
Wood, p.36
legend("topright", legend = levels(water$location), pch = c(1,2), bty = "n")
S-Plus であれば，marks 。
* 重複を取る [#l940a327]
例:/home/ishida/research/corpus/file.old/bak/file/writers19.R を参照
ここで unique は重複を取る。
legend(-4.3, 4.5, unique(writers.name),
fill=(as.integer(unique(writers.name))), text.width=1)
* 行列・配列の特定の次元を軸にあるベクトルを順に操作する、関数 sweep() [#x7d39075]
> x = matrix(1:16, 4,4)
> x
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 2 6 10 14
[3,] 3 7 11 15
[4,] 4 8 12 16
> y= 1:4 # sweep されるベクトル
> sweep(x,1,y) # x の各列から y を引き去る(既定動作)
[,1] [,2] [,3] [,4]
[1,] 0 4 8 12
[2,] 0 4 8 12
[3,] 0 4 8 12
[4,] 0 4 8 12
> sweep(x,2,y) # x の各行から y を引き去る(既定動作)
[,1] [,2] [,3] [,4]
[1,] 0 3 6 9
[2,] 1 4 7 10
[3,] 2 5 8 11
[4,] 3 6 9 12
> sweep(x,1,y, "*") # x の各列に y を掛ける
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 4 12 20 28
[3,] 9 21 33 45
[4,] 16 32 48 64
* updateの使い方 2005 /06 /04 [#l829ce4d]
例えば
writers.lda table(writers$writer,predict(writers.lda)$class)
と解析した結果と，und を取り除いた結果
writers.lda.und table(writers$writer,predict(writers.lda.und)$class)
を比較する。なお以下はエラーとなる
#writers.lda.step #anova(writers.lda , writers.lda.und)
* 判別分析における交差妥当化 2005/06/04 [#ub22b09e]
例えば
writers.lda #table(writers$writer,writers.lda$class) エラー
table(writers$writer,predict(writers.lda)$class)
#cross validation
writers.lda.cv CV=T)
table(writers$writer,writers.lda.cv$class)
#table(writers$writer,predict(writers.lda.cv)$class) エラー
* グラフに直線を引く 2005 / 06 / 14 [#j93a5c1d]
plot(x ~ y)
kaiseki.lm(x ~ y)
abline(kaiseki.lm[1,1], kaiseki.lm[3,1])
として切片，傾きを指定しても，また
abline(kaiseli.lm)
とするのも可能。
* 分散分析の結果に summary.lm を実行 2005/06/15 [#ff0ccd41]
Crawely p.163
分散分析で，平均値に差があるかを確認し，さらに，要因レベルの効果差を見るには，
summary.lm(aov(ozone ~ garden))
などを実行する。
* quartile と quantile percentile の使い分け 2005/06/15 [#odf1c97d]
Dalgaard p.58
//0.25, 0.5, 0.75 点の三つを three quartile といい，
//これに min, max を含めてそれぞれを quantile という。
Crawely p.197
//50% - 75% quartile と範囲を指定
//75% quantile などと分位点を指定
//
//quantile は，四分位点を一般化した p分位点，すなわち，100pパーセント点(0
* Boxplot での notch＝Ｔ オプション 2005/06/15 [#aee66ae1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.167
//ノッチは，メディアンの中心線から，上下に同じ長さの裾を伸ばすが，それぞれの限界値で折り曲がることによって，幅を示す。
* treatment.control の見方 2005/06/16 [#i78e2ed7]
[[Crawley Statistics Using R:http://www.amazon.co.jp/Statistics-Introduction-Michael-J-Crawley/dp/0470022981/]] p.174, [[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.287
treatment.control では，各要因の各レベルの効果を，
その要因の切片(基軸)であるレベルと比較して判断する。
したがって，基軸以外のレベルの相互の有意差は表示されていない。
基軸以外のレベル相互の比較
要因Aと他要因 BCD全体との差，BCD 内の差を一度に検定 2006 03 20
[[John Fox:http://socserv.mcmaster.ca/jfox/]], p.142 - 143.
他に[[ここ>#ContrastsMatrix]]
* loglin2 2*2 分割表の Likelihood ratio を求める [#u3b5389d]
&aname(loglin2);
始めに ここを参照
こことここも
またここも
meyer p.132 の例より
lessWords moreWords
USA 57 2
Phil 71 12
NZ 66 16
GB 23 0
meyerP.132 dimnames(meyerP.132)
model0 summary(model0)
# anova(model0) # これは用意されていない
model1 summary(model1)
model1$lrt
#P値は
1- pchisq(model1$lrt, model1$df)
1- pchisq(model1$pearson, model1$df)
* コレスポンデンス分析 [#q882e6d5]
R multiv パッケージのコレスポンデンス分析で出力される
行座標 $rproj, 列座標 $cproj はデフォルトでは，求められた x,y に単相関係数を乗じている。(言い方を変えると，固有値を平方根を乗じている)
参照 高橋信「Excelによるコレスポンデンス分析」, p.127
/home/research/statistics/correspondence.R を参照
#高橋信「Excelによるコレスポンデンス分析」
Z 13,8,15,16,
18,11,14,8), ncol=4, byrow=T)
#行ラベルとして使用する第1列をrowvarに代入
Z.rowvar<- c("M","H","C")
#rowvarの値を行ラベルとして読み込む
Z.colvar
rownames(Z) colnames(Z)
> Z
A B C D
M 10 19 13 5
H 13 8 15 16
C 18 11 14 8
Z.diagC #> Z.diagC
# [,1] [,2] [,3] [,4]
#[1,] 41 0 0 0
#[2,] 0 38 0 0
#[3,] 0 0 42 0
#[4,] 0 0 0 29
#>
Z.diagR #> Z.diagR
# [,1] [,2] [,3]
#[1,] 47 0 0
#[2,] 0 52 0
#[3,] 0 0 51
#>
Z.diagR.sqrt # [,1] [,2] [,3]
#[1,] 6.855655 0.000000 0.000000
#[2,] 0.000000 7.211103 0.000000
#[3,] 0.000000 0.000000 7.141428
#>
Z.solve Z.eigen
#$values
#[1] 1.00000000 0.07629421 0.01828369
#$vectors
# [,1] [,2] [,3]
#[1,] -0.5597619 0.74374060 0.3653992
#[2,] -0.5887841 -0.66725782 0.4561802
#[3,] -0.5830952 -0.04021102 -0.8114081
#>
# [,1] [,2] [,3]
# [1,] 0.3579767 0.2947642 0.3186919
# [2,] 0.2947642 0.3844402 0.3385965
# [3,] 0.3186919 0.3385965 0.3521610
# 固有値 １ と，対応する固有ベクトルは仮定
# 平均０，分散１に反する
# p.121 # x, y が求まる，以下はxのみ
Z.kai.X [,1] [,2] [,3]
[1,] -1 1.32867325 -0.6527762
[2,] -1 -1.13328105 -0.7747835
[3,] -1 -0.06896133 1.3915534
# 以下は Y のみ
Z.kai.Y # [,1]
# [1,] -0.23728726
# [2,] 1.46910930
# [3,] -0.05964337
# [4,] -1.50318462
Z.kai.Y # [,1]
# [1,] 1.5238380
# [2,] -0.6410599
# [3,] -0.1102451
# [4,] -1.1547168
# 固有値の平方根 は単相関係数
sqrt(Z.eigen$values)
# [1] 1.0000000 0.2762141 0.1352172
# p.127
Z.kai[,2] * sqrt(Z.eigen$values)[2]
# [1] 0.36699824 -0.31302816 -0.01904809
Z.kai[,3] * sqrt(Z.eigen$values)[3]
# [1] -0.08826657 -0.10476405 0.18816195
#(solve(Z.diagR.sqrt, sqrt(150) * Z.eigen$vectors)) %o% sqrt(Z.eigen$values)
> library(multiv)
> Z.results> Z.results
$evals
[1] 0.07629421 0.01828369
$rproj
Factor1 Factor2
[1,] 0.36699824 -0.08826657
[2,] -0.31302816 -0.10476405
[3,] -0.01904809 0.18816195
$cproj
Factor1 Factor2
[1,] -0.06554208 0.20604911
[2,] 0.40578865 -0.08668233
[3,] -0.01647434 -0.01490704
[4,] -0.41520073 -0.15613757
$rcntr
Factor1 Factor2
[1,] 0.553150080 0.1335166
[2,] 0.445232993 0.2081003
[3,] 0.001616926 0.6583831
$ccntr
* 図の保存 eps 形式 2005 06 24 [#o1d0f24b]
dev.copy2eps(file="log.eps")
* 対数尤度比をコーパスで求める場合 2005 07 01 [#j65238fa]
Dunning の式は，Oakes p.172 にまとめられている。必要となるのは a,b,c,d，と
当該コーパスにおける総語数 N (結局 a+b+c+d に等しい)。
a = あるキーワードにたいして設定したスパンの内部に，その共起語が現れた回数
b = キーワードにたいして，その共起語以外の語がスパンに現れた回数
c = キーワード以外の語にたいして，その共起語がスパンに現れた回数
d = キーワード以外の語にたいして，その共起語以外の語が現れた回数
* 東大「人文社会科学の統計学」，飽和モデル，独立モデル 2005 07 13 [#rc7b5e81]
始めに ここを参照
こことここも
またここも
&aname(social.mobility2);
データは東大「人文社会科学の統計学」 p.277
social.mobility2<-
array(c(
517, 244, 250,
152, 351, 328,
7, 13, 140,
441, 159, 272,
130, 227, 351,
24, 21, 268
),
dim = c(3,3, 2),
dimnames =
list(
father = c("white","blue","agri" ), # = c(1)
child = c( "white","blue","agri"), # = c(2)
year= c("1985", "1975") # = c(3)
)
)
class(social.mobility2 )<- "table"
ftable(social.mobility2)
library(MASS)
# 以下は，3次の効果を外した同時独立モデルだが，ここに有意差は無い。すなわち year の効果は無い。
# 3次の効果を含めた飽和モデルは，自由度ゼロなので必ず適合する。p. 271
social.mobility2.loglm3<- loglm(~father + child + year + father:child + father:year + child:year, data=social.mobility2)
anova(social.mobility2.loglm3, test="Chisq")
同じことだが，以下でもよい
social.mobility2.loglm3<- loglm(social.mobility2, list(c(1), c(2), c(1,2), c(2,3), c(1,3)))
1 -pchisq(kekka$lrt,kekka$df)#
1 -pchisq(kekka$pearson,kekka$df)#
# 東大「人文社会の統計学」 p.276
# log m = u... + ui.. + u.j. + u..k + uij. + ui.k + u.jk (9.33)
# ここで uij. ui.k u.jk がそれぞれ，father:child, father:year, child:year
# を表しているが
# uijk (p.276下から12行目), すなわち father:child:year
# は含まれていないモデル
# なお u... は全体平均 grandmean は暗黙のうちに含まれる
# ui.. 行平均 u.j. 列平均 u..k 軸平均は それぞれ要因 father, child, year
# を指定して含まれる
* 分散分析での母数モデル，変量モデル，混合モデル(豊田秀樹) 2005 07 27 [#h688e93f]
S-Plus での混合モデルの構築
[[Crawley:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.673 -
ちなみに
# p.678
lmedata3.model2 で。
random=~1|Genotype
は Genotype の intercepts をランダム要因と見なす指定, p.678
# 豊田秀樹 違いを見抜く統計学 p.120, p.127, toyoda.Bunsan.R を参照
p120<- data.frame(type = rep(c("sea","mount"), each = 8),
develop = rep(rep(c("old","new"), c(4,4)),2) ,
region = rep(c("izu","kuju","karui","yatsu"), each=4),
price = c(18.0, 15.0, 12.0, 20.7, 7.6, 12.3, 15.0, 15.2, 11.7,
11.6, 14.5, 15.4, 7.8, 7.5, 4.2, 6.4))
#p120<- data.frame(type = rep(c("sea","mount"), each = 8),
develop = rep(rep(c("old","new"), c(4,4)),2) ,
price = c(18.0, 15.0, 12.0, 20.7, 7.6, 12.3, 15.0, 15.2, 11.7,
11.6, 14.5, 15.4, 7.8, 7.5, 4.2, 6.4))
#> p120
# type develop region price
#1 sea old izu 18.0
#2 sea old izu 15.0
#3 sea old izu 12.0
#4 sea old izu 20.7
#5 sea new kuju 7.6
#6 sea new kuju 12.3
#7 sea new kuju 15.0
#8 sea new kuju 15.2
#9 mount old karui 11.7
#10 mount old karui 11.6
#11 mount old karui 14.5
#12 mount old karui 15.4
#13 mount new yatsu 7.8
#14 mount new yatsu 7.5
#15 mount new yatsu 4.2
#16 mount new yatsu 6.4
attach(p120)
交互作用を含んだモデル，交互作用棄却後のプーリングモデル
p120.aov<- aov(price ~ type * develop)
summary(p120.aov)
p.123と同じ出力# 豊田秀樹 違いを見抜く統計学 p.120, p.127
type develop の結果だけ見れば，プーリングした結果と同じ(p.123下)
Df Sum Sq Mean Sq F value Pr(>F)
type 1 84.181 84.181 10.1788 0.007769 **
develop 1 115.026 115.026 13.9084 0.002877 **
type:develop 1 8.556 8.556 1.0345 0.329170
Residuals 12 99.242 8.270
summary.lm(p.120.aov)
> 84.181 / 8.270# AB の平均平方ではなく，誤差の平均平方で割っている, p.108
[1] 10.17908
> 115.026 / 8.270
[1] 13.90883
> 1 - pf(10.1788, 1, 12)
[1] 0.007769413
交互作用を変量としたモデル，反復測定
p121.aov = aov(price ~ type + develop + Error(type/develop))
summary(p121.aov)
p.121 と同じ出力
母数因子だが，一つの標本から繰り返し測定した
すなわち，交互作用は変量モデル
Error: type:develop
Df Sum Sq Mean Sq F value Pr(>F)
type 1 84.181 84.181 9.8392 0.1965
develop 1 115.026 115.026 13.4444 0.1695
Residuals 1 8.556 8.556
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 12 99.242 8.270
ネストを仮定したモデル
p127.aov summary(p127.aov)
p.127 と同じ出力
すなわち，地域はタイプにネストしている
Error: region:type
Df Sum Sq Mean Sq F value Pr(>F)
type 1 84.181 84.181 1.3624 0.3635
Residuals 2 123.581 61.791
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 12 99.242 8.270
ここで
> 61.791 / 8.270
> 1 - pf(7.471705 , 2, 12)
[1] 0.007804986
三段のネストを仮定したモデル
これは上の２段ネストと同じ結果になっている
p.157 とは異なり，残差が二つに分離されていないx
p.157.aov > summary( p.157.aov )
Error: region:type
Df Sum Sq Mean Sq F value Pr(>F)
type 1 84.181 84.181 1.3624 0.3635
Residuals 2 123.581 61.791
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 12 99.242 8.270
また crawley, S-Plus, p.213下を参照。
p120.aov # = aov(price ~ type/develop)
# = aov(price ~ type/develop + Error( type/develop))
summary(p120.aov)
#
# p.127 と同じ出力
#
# Df Sum Sq Mean Sq F value Pr(>F)
#type 1 84.181 84.181 10.1788 0.007769 **
#type:region 2 123.581 61.791 7.4715 0.007806 **
#Residuals 12 99.242 8.270
#anova(p120.aov) #Error
#
# Df Sum Sq Mean Sq F value Pr(>F)
#type 1 84.181 84.181 10.1788 0.007769 **
#develop 1 115.026 115.026 13.9084 0.002877 **
#type:develop 1 8.556 8.556 1.0345 0.329170
#Residuals 12 99.242 8.270
#summary.lm(p120.aov)
#> 84.181 / 8.270 # AB の平均平方ではなく，誤差の平均平方で割っている, p.108
#[1] 10.17908
#> 115.026 / 8.270
#[1] 13.90883
#> 1 - pf(10.1788, 1, 12)
#[1] 0.007769413
# エラ−モデル 1
p120.aov summary(p120.aov)
#
# p.121 と同じ出力
#
#Error: type:develop
# Df Sum Sq Mean Sq F value Pr(>F)
#type 1 84.181 84.181 9.8392 0.1965
#develop 1 115.026 115.026 13.4444 0.1695
#Residuals 1 8.556 8.556
#Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
#Residuals 12 99.242 8.270
# エラ−モデル 2
p120.aov summary(p120.aov)
#
# p.127 と同じ出力
#
#Error: region:type
# Df Sum Sq Mean Sq F value Pr(>F)
#type 1 84.181 84.181 1.3624 0.3635
#Residuals 2 123.581 61.791
#Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
#Residuals 12 99.242 8.270
# ネストモデル
p120.aov summary(p120.aov)
#
# p.127と同じ出力
#
# Df Sum Sq Mean Sq F value Pr(>F)
#type 1 84.181 84.181 10.1788 0.007769 **
#type:region 2 123.581 61.791 7.4715 0.007806 **
#Residuals 12 99.242 8.270
#anova(p120.aov) #Error
* 対数尤度を求める関数 2005 08 01 [#o1feec7a]
logLik, (crawley, S-Plus, p.121)
* リストへのアクセス 2005 08 03 [#r8e750c9]
test
> test
Error: factory
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 6 260.70 43.45
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
treats 4 8.324 2.081 0.4637 0.7617
Residuals 24 107.704 4.488
以上のリストで，Error の F value にアクセスするには
summary(aov(rate ~ treats + Error(factory)))[[2]][[1]]$"F value"[1]
と実行する。
Crawely.S.plus.R を参照
>
d for(j in 1:1000){
treats for(i in 1:7) treats treats d[j] }
* グラフィックス表示を一枚ずつ見る 2005 08 05 [#d0392771]
par(ask=TRUE)
* 共分散分析, Ancova モデル式におけるアスタリスク 2005 08 08, 2006 01 11 [#hc1a7d93]
difference in the slop は 交互作用の項を見る。
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.295
covariate については p.293
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.216, p.341
w が sex カテゴリーの変数，x が連続量の変数の場合
y ~ w * x
は，w の二つのカテゴリそれぞれに回帰式(intecept と slope) regression slope を推定しろという命令となる。
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.126
// S-plus
* [S] EES の評価途中にスペースがある [#yb52040f]
式に改行やスペースが挟間っていると，そこで処理が止まってしまう。
* [S] ESS で操作中で画像の選択がある [#t7d79e80]
例えば plot を含めて，式を一括評価する時
par(ask=T)
と設定した後だと，yes か no の答えを得ることができないので，そこで止まってします。
* データの並び替え [#q9c7bb32]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.417のデータならば
> fishery
xv yv
1 13.9885004 4.56836585
2 89.7028342 0.29710071
3 9.3151971 3.90699292
4 9.2132163 4.38109971
5 34.7017719 1.85936367
のようなデータフレームならば
fishery[order(fishery$xv),]
とすれば，xv の昇順にデータが整列する。
* [S] 内積とノルム 2005 08 23 [#f3b0090e]
a b
内積は
t(a) %*% b
としても
sum(a * b)
としても良い。
ノルムは
vecnorm(a), vecnorm(b) で求めるか
mynorm if(p=="inf") return(max(abs(x)))
sum(abs(x)^p)^(1/p))
}
で求める。
* adj の設定 2005 08 25 [#ub65cb39]
research/tex/p2005/2005/dunne.R より
text(E.normalized.B.dunne2[1, E.normalized.B.dunne2[2,] > 0.0],
E.normalized.B.dunne2[2, E.normalized.B.dunne2[2,] > 0.0],
adj=c(1.5,-.5), labels=( colnames(
E.normalized.B.dunne2[, E.normalized.B.dunne2[2,]
adj=c(X,Y) X を正にすると，点は左に，Yを正にすると点は下に移動
* 負の2項分布 2005 08 30 [#w1336322]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.483
例えば
成功率 0.6, 成功までの試行数が 5 として(すなわち６回目に成功)，
100回の試行を行なう。
count
table(count)
0 1 2 3 4 5 6 7 8 9
8 13 20 19 18 7 9 3 2 1
この場合 失敗無しで５回の成功が出た試行が 8 回
失敗 1 で５回の成功が出た試行が 13 回
あるのを意味する。
# 2006 07 24 追加
工学のためのデータサイエンス, p.49 - p.50
R の関数 _nbinom の引数 size, prob は
平均 mean が 'size * (1 - prob)/prob' を満たし，また
確率 prob が 'scale/(1 + scale)' を満たすように取られる
d p(0, 1612), rep(1,164), rep(2,71), rep(3,47), rep(4,28), rep(5,17),
rep(6,12), rep(7,12), rep(8,5), rep(9,7), rep(10,6), rep(11,3), rep(12,3),
rep(13, 13))
mean(d) # = 0.6005
var(d) # = 3.262531
prob = 0.1155 * (1 - 0.1553) / 0.1553# = 0.6282218 = size *(1 - prob)/prob
scale = 0.1838523 / ( 1 + 0.1838523)# = 0.1553 = scale/ (scale + 0.1553)
* lm, glm オブジェクトの summary における t 検定,z 検定 2005 09 02 [#t54cd3ab]
S-Plus では t 値が，R では z 値がデフォルト
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.534
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] S-plus.R, p.534
また S-Plus では，lm の出力は Pr(|t|) を出力するが，glm ではデフォルトでは出力されない。
* rgl.bbox の引数 2005 09 02 [#s091d17e]
rgl.bbox(color="#999999", emission="#FFA500",specular="#FFF0F5",
shinines=5,alpha=0.8)
color="#999999" は座標の数値ラベル色
emission="#FFA500" はボックスの壁の色
specular="#FFF0F5" は明るさ，数値が大きいほど明るい
shinines=5 は反射の具合
* 関数グラフの書き方 2005 09 19 [#q2e33709]
x y plot(x,y,type="l")
同じことだが
curve((x^2 - 6*x + 10), 0,6)
* 2項分布を Poisson 分布で近似する 2005 09 28 [#h5577017]
例えば2項分布で，0.00001 で成功する賭けを 100000 回行なった時，
ちょうど２０回成功する確率は
pbinom(20, 100000, .0001) - pbinom(19, 100000, .0001)
[1] 0.001865335
これをポアソン分布で近似するには
ppois(20, 10) - ppois(19, 10)
[1] 0.001866081
* データフレームのコピー 2005 10 08 [#s5f3e4ac]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.63 - 64
> search()
[1] ".GlobalEnv" "package:car" "package:methods"
[4] "package:stats" "package:graphics" "package:grDevices"
[7] "package:utils" "package:datasets" "Autoloads"
[10] "package:base"
> cooperation
エラー：オブジェクト "cooperation" は存在しません
> attach(Guyer) # Guyerデータには変数 operationがある
> objects()
[1] "Guyer" "hump.data" "hump.model" "last.warning"
> cooperation
[1] 49 64 37 52 68 54 61 79 64 29 27 58 52 41 30 40 39 44 34 44
> cooperation > objects()新たに変数 cooperation が作成されている
[1] "Guyer" "cooperation" "hump.data" "hump.model" "last.warning"
> rm(cooperation)# 新たに作成された cooperation を消し
> cooperation[1] Guyer変数のcooperation[1]を変更したつもりでも
> objects()実際にはコピーができている
[1] "Guyer" "cooperation" "hump.data" "hump.model" "last.warning"
> rm(cooperation)
> objects()
[1] "Guyer" "hump.data" "hump.model" "last.warning"
* png などの画像への保存 2005 11 09 [#b2a9096e]
始めに
png(file="gazo.png") #として device を平き
polt(y ~ x) # として画像を作成し
dev.off() # として device を閉じる。
ここで始めて保存される。
* 標準誤差 2005 12 02 [#zdcdd38c]
標準誤差とは，平均値の推定値の範囲
標準偏差とは，データが平均を中心にどれだけばらついているかの目安
* Boxplot の notch の解釈 2005 12 23 [#mcf9c930]
Boxplot の箱は 75 - 25 %タイルを示しているが，あるデータのこの範囲が，他のデータの平均値と重なっていない場合，二つのmedian は有意に異なる。Tukey による
[[Crawley Statistics Using R:http://www.amazon.co.jp/Statistics-Introduction-Michael-J-Crawley/dp/0470022981/]] p.167
* [S] 関数内で式を表示させない 2005 12 24 [#e11ee55f]
関数内では，付置の式以外は表示される。
特に print 文では注意が必要である。
これを防ぐため，invisibleを使う。
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/statcomp/]] p.35
渋谷，柴田 p.34
* [S] S-plus + ESS 環境でのヘルプ 2005 12 24 [#d40dc98b]
キーワードにカーソルを合わせ，C-c C-v で表示させる。
?キーワードでは表示されない
// <<- はグローバル変数への代入
//
// 舟尾 p.183