トップ
新規
一覧
単語検索
最終更新
ヘルプ
ログイン
アールメカブ
R_fromOldHtml1
をテンプレートにして作成
開始行:
個人メモを整理中です．不適切な記述が多々あるかと思います...
#contents
* [[ ]] の使い方 [#z6109ecf]
例えば以下のようなデータで ([[Maindonald 旧版:http://wwwm...
> 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)
>
ちなみに，外れ値はレンジで指定するが，デフォルトでは rang...
またrug() で，y軸右(side=4)ないし左(side=2)に，プロットに...
ラグプロットについては[[Rjp Wiki:http://www.okada.jp.org/...
なお boxplot と split, rug を併用した例は [[Maindonald 旧...
* リストの先頭要素を取る [#r796e7ac]
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()は分位点関数なので，正規...
> 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-boo...
X11() を実行。サイズを指定するには X11(width=7.5,height=4)
現在アクティヴなグラフィックデバイスを閉じるには
dev.off()
特定のデバイス, 例えば ３ 番を閉じるには
dev.off(3)
デバイスの変更, 例えば ２ 番に切替えるなら
dev.set(2)
* 特定の行の取り出し [#l372d554]
例えば [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~john...
> 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
1 A Sprayed Slugs Fenced Unlimed Control
...
のようなデータフレームがあるとして，Block が A の項目は
splitplot[splitplot$Block=="A",]
として得られる。
~/mysplus/crawley.s , p.352
次の例は，term.cos で，行名が Berg である行の各列から，0....
[berg >= 0.999]
# あるマトリックスの，ある名前の行の (1列目の) 数値を見る
stif.keller.without.svd.Doc3
[rownames(stif.keller.without.svd.Doc3)=="Gouverneur"...
* 図における曲線，座標設定などの詳細 [#j0fca1af]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
なお plot で引数に type="n"を指定すると，点は刻印されない...
points() を実行する。
二つのカテゴリについて，変数 Y と 変数 X の対応をプロット...
* aggregate の使い方 [#cbfa9ae1]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
例えば
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 # このままだと元のデータフレムの情報として，Specie...
levels(iris.2$Species)
#[1] "setosa" "versicolor" "virginica"
そこで三つ目のレベルを欠損値とする。これによりレベルは二...
levels(iris.2$Species)[3]
iris.nnet iris.nnet.pre table(iris.2$Species, iris.nnet...
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 の固有値を求める。このと...
* 自己組織化マップ [#pffef85e]
[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?cmd=re...
#
library(som)
library(class)
ex gr#ex[1,]
testplot(test)
# 六角形配置にしたがって円を配置する (下記の 12*8は xdim,...
symbols(test$grid$pts[,1],test$grid$pts[,2],
circles=rep(0.5,8*12),inches=F,add=T)
# knn1を使って ex 中のデータが SOMで割り当てた test$code ...
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(...
と計算される。
* 作業ディレクトリにあるすべてのオブジェクトの削除 [#u8b...
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-boo...
> 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]
[[旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] C...
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-boo...
* set.seed [#peb6a616]
set.seed(x) は，R起動には，x をランダムな値に設定して生成...
したがって，同じ x の値で実行すれば，同じ乱数のセットが得...
結果として，同じ結果が得られる。[[Maindonald 旧版:http://...
* split [#s1d4876e]
第一引数で指定したベクトルを，第二引数で指定したカテゴリ...
Maindonald p. 134, p.213
* text関数の使い方 [#a0c82837]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
plot(weight.split$hb ~ volume.split$hb,
pch=16, xlim=range(volume), ylim = range(weig...
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 がマイナスなら，左に，また...
* 回帰分析で，変数間の相関係数を出すオプション [#f7ea1c56]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
allbacks0.lm summary(allbacks.lm, corr=T)
また個々のサンプルを抜いた場合の影響度を調べるには [[Main...
round(coef(allbacks.lm0), 2)
round(lm.influence(allbacks.lm0)$coef, 2)
* 変数間の線型性を求める alias [#j5f121ed]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
* 分散分析，回帰分析で base ラインを変更する relevel [#s...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
* option の使い方 [#e595e701]
Maindonals p.178
* 散布図の外にプロットを打つ points(, ..., xpd =T) [#cf9...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
* nlme マルチレベルでの固定効果と，ランダム効果 [#w367d...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
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]]と[[こ...
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 ...
飽和モデルの結果をみる
# 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]]も．[[...
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 + Fo...
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
>
##### よって最適なモデルは
# ~Length + Form + Type + Length:Form + Form:Type
#####
別データによる分析
それぞれの効果の推定値には，ダミーコーディングによる場合...
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","...
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....
(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...
の表3.2「ロケットの発射試験」
この例題の飽和モデルによる推定値は
http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03...
の表3.11にダミーコーディング、表3.12にANOVAコーディング
なお SPSS で出力される Partial association は「標準化係...
以上は[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?...
終了行:
個人メモを整理中です．不適切な記述が多々あるかと思います...
#contents
* [[ ]] の使い方 [#z6109ecf]
例えば以下のようなデータで ([[Maindonald 旧版:http://wwwm...
> 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)
>
ちなみに，外れ値はレンジで指定するが，デフォルトでは rang...
またrug() で，y軸右(side=4)ないし左(side=2)に，プロットに...
ラグプロットについては[[Rjp Wiki:http://www.okada.jp.org/...
なお boxplot と split, rug を併用した例は [[Maindonald 旧...
* リストの先頭要素を取る [#r796e7ac]
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()は分位点関数なので，正規...
> 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-boo...
X11() を実行。サイズを指定するには X11(width=7.5,height=4)
現在アクティヴなグラフィックデバイスを閉じるには
dev.off()
特定のデバイス, 例えば ３ 番を閉じるには
dev.off(3)
デバイスの変更, 例えば ２ 番に切替えるなら
dev.set(2)
* 特定の行の取り出し [#l372d554]
例えば [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~john...
> 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
1 A Sprayed Slugs Fenced Unlimed Control
...
のようなデータフレームがあるとして，Block が A の項目は
splitplot[splitplot$Block=="A",]
として得られる。
~/mysplus/crawley.s , p.352
次の例は，term.cos で，行名が Berg である行の各列から，0....
[berg >= 0.999]
# あるマトリックスの，ある名前の行の (1列目の) 数値を見る
stif.keller.without.svd.Doc3
[rownames(stif.keller.without.svd.Doc3)=="Gouverneur"...
* 図における曲線，座標設定などの詳細 [#j0fca1af]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
なお plot で引数に type="n"を指定すると，点は刻印されない...
points() を実行する。
二つのカテゴリについて，変数 Y と 変数 X の対応をプロット...
* aggregate の使い方 [#cbfa9ae1]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
例えば
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 # このままだと元のデータフレムの情報として，Specie...
levels(iris.2$Species)
#[1] "setosa" "versicolor" "virginica"
そこで三つ目のレベルを欠損値とする。これによりレベルは二...
levels(iris.2$Species)[3]
iris.nnet iris.nnet.pre table(iris.2$Species, iris.nnet...
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 の固有値を求める。このと...
* 自己組織化マップ [#pffef85e]
[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?cmd=re...
#
library(som)
library(class)
ex gr#ex[1,]
testplot(test)
# 六角形配置にしたがって円を配置する (下記の 12*8は xdim,...
symbols(test$grid$pts[,1],test$grid$pts[,2],
circles=rep(0.5,8*12),inches=F,add=T)
# knn1を使って ex 中のデータが SOMで割り当てた test$code ...
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(...
と計算される。
* 作業ディレクトリにあるすべてのオブジェクトの削除 [#u8b...
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-boo...
> 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]
[[旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] C...
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-boo...
* set.seed [#peb6a616]
set.seed(x) は，R起動には，x をランダムな値に設定して生成...
したがって，同じ x の値で実行すれば，同じ乱数のセットが得...
結果として，同じ結果が得られる。[[Maindonald 旧版:http://...
* split [#s1d4876e]
第一引数で指定したベクトルを，第二引数で指定したカテゴリ...
Maindonald p. 134, p.213
* text関数の使い方 [#a0c82837]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
plot(weight.split$hb ~ volume.split$hb,
pch=16, xlim=range(volume), ylim = range(weig...
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 がマイナスなら，左に，また...
* 回帰分析で，変数間の相関係数を出すオプション [#f7ea1c56]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
allbacks0.lm summary(allbacks.lm, corr=T)
また個々のサンプルを抜いた場合の影響度を調べるには [[Main...
round(coef(allbacks.lm0), 2)
round(lm.influence(allbacks.lm0)$coef, 2)
* 変数間の線型性を求める alias [#j5f121ed]
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
* 分散分析，回帰分析で base ラインを変更する relevel [#s...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
* option の使い方 [#e595e701]
Maindonals p.178
* 散布図の外にプロットを打つ points(, ..., xpd =T) [#cf9...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
* nlme マルチレベルでの固定効果と，ランダム効果 [#w367d...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
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]]と[[こ...
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 ...
飽和モデルの結果をみる
# 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]]も．[[...
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 + Fo...
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
>
##### よって最適なモデルは
# ~Length + Form + Type + Length:Form + Form:Type
#####
別データによる分析
それぞれの効果の推定値には，ダミーコーディングによる場合...
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","...
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....
(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...
の表3.2「ロケットの発射試験」
この例題の飽和モデルによる推定値は
http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03...
の表3.11にダミーコーディング、表3.12にANOVAコーディング
なお SPSS で出力される Partial association は「標準化係...
以上は[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?...
ページ名:
