トップ
新規
一覧
単語検索
最終更新
ヘルプ
ログイン
アールメカブ
R_fromOldHtml2
をテンプレートにして作成
開始行:
個人メモを整理中です．不適切な記述が多々あるかと思います...
#contents
* 互いに相関しているサンプルの分散の差の検定 2006 01 06 ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* 偏相関係数の計算方法 2006 01 06 [#w512ce2d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
lm で，対象としたい変数を取り除いたモデルを update で生成...
その変数の 偏 SSR (regression sum of squares) が求まる。
これを SST で割り，平方に開けば良い。
* SST (total sum of squares) の計算方法 2006 01 06 [#z45...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
lm(X ~ 1)
* Anova 表の見方，RSS 残差平方和 2006 01 06 [#xe0ae3d7]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
RSS は回帰によって説明されない変動，長谷川「多変量」 p.22
* [S] S-Plus で locale のエラー [#bb3c3fa8]
Warning: Locale not supported for XmbTextListToTextProp...
Warning: Cannot convert XmString to compound text
http://s.isac.co.jp/FAQ/faq-0032.html を参照
http://S.isac.co.jp/docs-4.5/machdep_3.htmlを参照
その他，myTurboAMD.html の S-plus インストールを参照
* [S] detach rm の順番 2006 01 09 [#qd4b3d6a]
まず detach し，次に rm する。逆に行なうと，
Problem: Object "factories" not found
Use traceback() to see the call stack
とエラー表示
* ネストの式 2006 01 11 [#b93e19cd]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
y ~ A/B
y ~ A * A:B
y ~ A + B %in % A
に等しい。
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.135 - 136
複数の要因からなるデータの nest は，ある要因の水準ごとに
別の要因の傾きを算出する。そして
A) すべての切片が 0 か
B) ある要因のある水準に対する別の要因の傾きがすべて 0 か
を検定する。
ネストモデルについては三中先生のサイト
http://cse.niaes.affrc.go.jp/minaka/R/NestedANOVA.html
を参照
* 切片を除く．その解釈 2006 01 11 [#t4e8647c]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p. 132 下, ...
y ~ x - 1
は説明変数が連続量の場合，切片を 0 にし，全体平均と 全体 ...
すべての説明変数がカテゴリデータの場合，各レベルの平均を...
( -1 が無ければ，全平均，また平均からの差になる)。
* scale location plot 2006 01 14 [#n7eb0bef]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
回帰分析の結果を plot の引数に与えた場合の第三のプロット...
分散が定数であるかを確かめるのに有効。
* カテゴリベクトルの順序づけ ordered 2006 01 14 [#f2988f...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.130
* 変数削減の方法 anova と summary の違い 2006 01 18 [#ad...
Dalgaard, S-Plus, p.181, p.155
summary() で表示された P 値は，表示順に関係なく削除可能。
これに対して anova() の場合は，表示順序が重要。
また高次の相互作用をチェックする場合，ある単独ファクター...
それを含む高次作用が有意の場合は，その単独ファクターは残...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* 適合値と分散の対応 2006 01 18 [#f519a964]
要因A と要因Bの各レベル語とに組合せにおける実験回数が少な...
多少，適合値に対して分散が減少する傾向があっても問題ない。
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* [S] AIC の計算方法 2006 01 27 [#h08b644e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
mode.1<- lm(growth ~ 1)
AIC(model)
* NULL model はパラメータを grand mean 以外推測しない 20...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
model<- aov(y ~ 1)
model<- lm(y ~ 1)
* Helmert contrasts vs. Treatment contrasts 2006 01 30 [...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
y = u + b1 x1 + b2 x2 + b3 x3 + b4 x4
というモデルが想定される場合，五つのパラメータを想定する...
treatment mean は
b1 を 0 と見なし，x1 の平均を u とし，他の b は，u との差...
これに対して Helmert Contrasts は，
全体平均を最初のパラメータとし
最初と二つ目の平均の平均と，最初の平均との差を第二のパラ...
最初と二つ目の平均の平均と，最初と二つ目さらに三つ目の平...
最初と二つ目さらに三つ目の平均の平均と，全体平均の差を第...
またコントラストを手作業で設定するのは [[John Fox:http://...
他にもここを参照
* 繰り返しのない二元配置分散分析 2006 01 31 [#af9d9e29]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Block, Plot, Split-plot の違い 2006 02 01 [#w63cb6bb]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* tapply と list の組み合わせ 2006 02 01 [#fcee3a67]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
attach(splityield)
names(splityiled)
> names(splityield)
[1] "yield" "block" "irrigation" "density"
"fertilizer"
だとして
tapply(yield, list(block, irrigation), sum)
等のように使う。
list については他に crawley, S-Plus, p.406, p.414
* Interaction が有意の時，main effect は削除できない 200...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Fixed Effects はそのファクターレベルに情報がある 2006 ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* [S] interaction.plot の使い方 2006 02 02 [#l37c66ab]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* [S] model.table の使い方 2006 02 02 [#s2e21fb1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Block と Plot, Split-Plot 2006 02 03 [#qc63bace]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
Block (= Plot) をさらに半分にして Split-Plot と呼ぶ。
* Treatment structure, Error structure 2006 02 03 [#efa1...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Predict と estimate の違い 2006 02 03 [#u93595c1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
また fitted と predict の違いは [[Faraway:http://www.amaz...
* Corrected, uncorrected Sum os Squares 2006 02 03 [#g33...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
ある変数の値 に その変数の値 をかけた値の総和 uncorrected...
から
その変数の総和の二乗をその変数の数で割ったものを引くと co...
Corrected, uncorrected Sum os Squares の違いは，前者は変...
* 変数変換について 2006 02 07 [#a1daf86a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* predict の使い方 2006 02 06 [#ra24fe2d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* sample with replacement 2006 02 08 [#v2aabfa4]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
sample with replacement, 取り出した石をもとに戻す
sample without replacement もとに戻さない Hypergeometric
* ヒストグラムで整数値をバーの中心にする 2006 02 09 [#d6...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
break points を 0.5 間隔にする。
hist(mass, breaks = -0.5:16.5)
* Deviance (乖離度) について 2006 02 11 [#x276d817]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
単独モデルの場合
total deviance と residual deviance は
それぞれ satured model との差を表している。
複数モデル比較の場合
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
「パラメータの少ない方のモデルがフィットしている」が帰...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
regression 等ではモデルの当てはめ具合を SSE で測ったが，
GLM では，binomial な proportion データについて deviance ...
また p.526
Deviance is equivalent to sums of squares in linear mod...
p.539 には
total deviance と residual deviance の違い。
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
に単独のモデルについての deviance と，
二つの入れ子になったモデル(どちらかがパラメータが多い)を...
* GLM 小サンプルの weight を小さくする 2006 02 14 [#n8ef...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
fail <- total - success
y <- cbind(success, fail)
model1 <- glm(y ~ explanatory,binomial)
* カテゴリ変数のレベル変更 (アルファベット順以外) 2006 0...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
ordered(target, levels = c("E","D","C","B","A"))
* binomial, exponential, poisson エラーについて(overdisp...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
Everit & Hothorn p.103 - 105.
dispersion parameter は 1 にセットされる。
なぜなら binomial, exponential, poisson では分散は，デー...
residual deviance が residual degrees of freedom よりも大...
overdispersion
この時，モデルの比較はデビアンスではなく，F検定量で行う (...
dispersion パラメータの推定方法は Faraway, p.47, p.60
* logistics analysis of covariance with proportion data ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
例えば三つの薬品のそれぞれについて，分量を変えながら，効...
* wilcox.testが有効ではないケース 2006 02 17 [#q7dce0b3]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
0,1 からのみなるようなデータの場合同順位が多くなる。
* G test，chisq とは別の方法による独立性の検定 2006 02 1...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
glm.crss<- glm(c(22,61,69) ~1, family=poisson)
1- pchisq( glm.crss[13]$deviance, glm.crss[7]$df.residu...
あるいは
glm.crss<- glm( c(15,10,10,15) ~ 1, poisson)
summary( glm.crss)
1 - pchisq(glm.crss$deviance, glm.crss$df.residual)
[1] 0.5695988]
これは
pchisq(glm.crss$deviance, glm.crss$df.residual,
lower = FALSE)
と同じこと
* nuisance 変数 2006 02 18 [#gb4dd871]
crawley, p.551
* 高次の交互作用がある場合，それより低次の要因は省けない...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* モデル比較を多数くり返す場合の p 値 2006 02 18 [#n9a36...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
"With this many comparisons, I would always work at
p = 0.01 for significance, rather than p = 0.05 "
* 二つの要因に有意差があると出力されても，それは Interce...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
summary(p.580.model)
Call: glm(formula = cbind(relief, patients - relief)
~ treat, family = binomial)
Coefficients:
Value Std. Error t value
(Intercept) -1.067841 0.2471428 -4.320744
treatB 1.959839 0.3427433 5.718094
treatC 2.468734 0.3666004 6.734128
# これは Intercept と差
* GAM (Generalised Additive Models 2006 02 24 [#x27432da]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
y ~ s(x) + s(y) + s(z)
y ~ s + lo(w) + z
s は s smoother , lo は lo smoother を表す。
* Trend がある 2006 02 28 [#v681bb1e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
回帰分析をして，傾きに有意差があれば，そこにトレンドが認...
* プロット描画領域にテキストを加える 2006 03 02 [#wa2f9c...
text 関数を使う。例えば
text(grexp.obs$N[nrow(grexp.obs) ] ,
grexp.obs$K[nrow(grexp.obs) ], labels = "grexp")
これはデータフレーム grexp.obs の N 列の最大値を X 軸に，...
baayen.R を参照。
* スケールの違うプロットを重ねる 2006 03 03 [#p0db652a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
par(mar = c(5,4,4,5) + 0.1) としてマージンを広げ
tsplot(cinnabar, ylab = "cinnabar") 最初のグラフを描き
par(new = T, xaxs = "d") として新しいプロットが先のプロッ...
* 空のリストを作成する方法 2006 03 03 [#n3bb2f43]
Ligges p. 103
X <- vector(mode = "list", length=10)
* リストに対する演算 2006 03 03 [#w35e23a7]
X <- vector(mode = "list", length=5)
# 5個の要素からなるリストを作成し
lapply(lapply(compare.K, "[", , 2 )
lapply(lapply(compare.K, "[", , "N" )
# リストの要素(データフレーム)のすべてから，
# すべての行の 2 列目(あるいはラベル名が N の列)を取り...
lapply(lapply(compare.K, "[", , 2 ) , "max")
# リストから上記の条件で抽出した各要素について，
# その最大値を取り出す
max(as.numeric (lapply(lapply(compare.K, "[", , 1 ) ,
"max") ))
# 最終的に，リスト全要素の 2 列目の最大値を取り出す
* sapply の使い方 2006 03 04 [#k8e3c932]
Z.list の各要素はベクトルだとする。
unlist(Z.list) でベクトルが返される。
ちなみに sapply(Z.list, "[",1)
は最初の要素をベクトルとして取り出す。
sapply(Z.list, "[") はリストの各要素(ベクトル)を列とした...
Z.data<- data.frame(sapply(Z.list, "["))
とすればデータフレームに変更可能
* paste, parse, expression, eval によるデータフレームの...
x <- "col"
y <- 1:5
z <- "th = 0"
xyz <- paste(x,y,z,sep="")
# xyz
# [1] "col1th = 0" "col2th = 0" "col3th = 0"
"col4th = 0" "col5th = 0"
xyz <- paste(xyz, collapse=",")
# xyz
# [1] "col1th = 0,col2th = 0,col3th = 0,
col4th = 0,col5th = 0"
d <- "data.frame("
d1 <- ")"
d2 <- paste(d,xyz,d1)
# d2
# [1] "data.frame( col1th = 0,col2th = 0,col3th = 0,
col4th = 0,col5th = 0 )"
parse(text =d2)
# parse(text =d2)
# expression(data.frame(col1th = 0, col2th = 0,
col3th = 0, col4th = 0,
# col5th = 0))
Z <- eval(parse(text =d2) )
# Z <- eval(parse(text =d2) )
# Z
# col1th col2th col3th col4th col5th
# 1 0 0 0 0 0
* 文字列の処理，文字列を式に変更 2006 03 04 [#s9534d8b]
Ligges, p.54, p.91
例えばデータフレームの列名をベクトルに変換するには
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
pig<- read.table(
"/S/crawley/data/pig.txt",header=T)
attach( pig )
names( pig )
deparse( names( pig ))
# "c(\"Pig\", \"t1\", \"t2\", \"t3\", \"t4\", \"t5\",
\"t6\", \"t7\", \"t8\", \"t9\")"
parse(text = deparse( names( pig )))
#expression(c("Pig", "t1", "t2", "t3", "t4", "t5",
"t6", "t7", "t8", "t9"))
eval(parse(text = deparse( names( pig ))))
#[1] "Pig" "t1" "t2" "t3" "t4" "t5"
"t6" "t7" "t8" "t9"
parse(text = eval(parse(text = deparse( names( pig )))))
#expression(Pig, t1, t2, t3, t4, t5, t6, t7, t8, t9)
よって次のように行なう
pig.colnames < - eval(parse(text = deparse(names(pig))))
* デ−タフレ−ムやリストをベクトルに変換 2006 03 06 [#n4d2...
x1<- data.frame(y1 = c(1,2,3), y2 = c(11,22,33))
x2<- data.frame(y1 = c(4,5,6), y2 = c(44,55,66))
x3<- data.frame(y1 = c(7,8,9), y2 = c(77,88,99))
X123<- list(x1,x2,x3)
リスト X123 の各要素はベクトルだとする。
unlist(X123) でベクトルが返される。
# データフレームの各列をつなげる
unlist(x1)
あるいは
attach(x1)
z<- NULL
for(i in 1:2){
z<- c(z, get(paste("y", i ,sep="")))
}
# 列数が少ないなら
c(y1,y2) #でも OK
* rep()関数の使い方 2006 03 06 [#w4f9f6e7]
> rep(c("Carroll","Dickens", "Eliot", "Bronte",
"Stevenson", "Gaskell"), rep(3,6))
あるいは
> rep(c("Carroll","Dickens", "Eliot", "Bronte",
"Stevenson", "Gaskell"),each = 3)
* 文字列操作 2006 03 06 [#m580f36e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
Ligges p.54
例えば strtrim(labs,2) は文字列ベクトル labs のそれぞれの...
* アルファベットの略語 2006 03 07 [#w7fcf51b]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
abbreviate(str)
* 対数尤度比を算出 ，松原「統計学の考え方」 2006 03 10 [...
松原 「統計学の考え方」, p82 - p.84
> ftable(p.73)
death Y N
suspect victim
Wh wH 19 132
bL 0 9
Bl wH 11 52
bL 6 97
m2 iterations: deviation 8.881784e-16
$lrt
[1] 0
$pearson
[1] 0
$df
[1] 0
$margin
$margin[[1]]
[1] 1
$margin[[2]]
[1] 2
$margin[[3]]
[1] 1 2
loglin(matu, c(1,2), param =T)
2 iterations: deviation 0
$lrt
[1] 0.1658080
$pearson
[1] 0.1690821
$df
[1] 1
あるいは
library(MASS)
matu<- array(data=c(4,6,20,40),
dim=c(2,2), dimnames = list(Row =c("a1","a2"),
Col = c("b1", "b2")) )
loglm(~Row + Col, data = matu)
x<- array(data=c(15, 10, 10, 15), dim=c(2,2))
# = x<- matrix(c(4,6,20,40), nrow = 2)
kekka<- loglin(x, c(1,2), param =T)
1 -pchisq(kekka$lrt,kekka$df)# 対数尤度比 G による検定...
1 -pchisq(kekka$pearson,kekka$df)# カイ二乗検定による検...
* R での時系列分析 2006 03 10 [#f9c56cde]
R-wiki
R の基本パッケージ base, stats 中の時系列オブジェクトの...
時系列データに関する最も基本的な統計量は
自己共分散(auto-covariance)と
自己相関係数(auto-correlation)である。
R の関連する関数は acf(), pacf(), ccf() である。
関数 acf() は時系列オブジェクトの自己共分散と自己相関係...
既定でそれをプロットする。
関数 pacf() は 偏自己相関係数(partial autocorrelations)...
関数 ccf() は二つの一次元時系列間のクロス相関係数
(cross-correlation)と
クロス共分散(cross-cavarinace)を計算する。
acf(x, lag.max = NULL,
type = c("correlation", "covariance", "partial"),
plot = TRUE, na.action = na.fail,
demean = TRUE, ...)
pacf(x, lag.max = NULL, plot = TRUE,
na.action = na.fail, ...)
ccf(x, y, lag.max = NULL,
type = c("correlation", "covariance"),
plot = TRUE, na.action = na.fail, ...)
* 回帰分析の summary の解釈 2006 03 14 [#c2d6ebc2]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.134
summary(prestige.mod.3)
Call:
lm(formula = prestige ~ income + education + type
+ income:type + education:type)
Residuals:
Min 1Q Median 3Q Max
-13.462 -4.225 1.346 3.826 19.631
Coefficients:
Estimate Std. Error t value Pr(>|t...
(Intercept) 2.276e+00 7.057e+00 0.323 0.74...
income 3.522e-03 5.563e-04 6.332 9.62e-...
education 1.713e+00 9.572e-01 1.790 0.07...
typeprof 1.535e+01 1.372e+01 1.119 0.26...
typewc -3.354e+01 1.765e+01 -1.900 0.06...
income:typeprof -2.903e-03 5.989e-04 -4.847 5.28e-...
income:typewc -2.072e-03 8.940e-04 -2.318 0.02...
education:typeprof 1.388e+00 1.289e+00 1.077 0.28...
education:typewc 4.291e+00 1.757e+00 2.442 0.01...
例えば，ここでは，Intercept は type要因の bc 水準で，これ...
income, education は type要因の bc 水準の時の傾き
typeprf, typewc は切片との差
income:typeprof は，type要因の prof 水準の時の傾き
Dalgaard, p.178
Call:
lm(formula = log10(diameter) ~ log10(conc) * glucose,
data = hellung)
Residuals:
Min 1Q Median 3Q Max
-2.672e-02 -4.888e-03 5.598e-05 3.767e-03 1.761e-02
Coefficients:
Estimate Std. Error t value Pr(>|t...
(Intercept) 1.627926 0.033754 48.230 < 2e-...
log10(conc) -0.046716 0.006846 -6.823 1.51e-...
glucose 0.003418 0.023695 0.144 0.8...
log10(conc):glucose -0.006480 0.004821 -1.344 0.1...
まず最初の二つは glucose があるデータの切片と傾き，次の二...
二つの切片の差は 0.003418 , 傾きの差は -0.006480
差は二つとも有意でない(0.88, 0.185)
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
options(contrasts=c("contr.treatment","contr.poly"))
Call: lm(formula = Weight ~ Age * Genotype)
Residuals:
Min 1Q Median 3Q Max
-0.7787 -0.3986 -0.01139 0.4055 0.8445
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 7.5407 0.3853 19.5724 0.0000
Age 0.2931 0.1162 2.5230 0.0150
GenotypeCloneB 1.0401 0.5449 1.9090 0.0623
GenotypeCloneC -0.9487 0.5449 -1.7411 0.0881
GenotypeCloneD 0.5883 0.5449 1.0797 0.2857
GenotypeCloneE -0.9587 0.5449 -1.7596 0.0849
GenotypeCloneF 1.5693 0.5449 2.8802 0.0059
AgeGenotypeCloneB -0.0241 0.1643 -0.1468 0.8839
AgeGenotypeCloneC -0.0316 0.1643 -0.1926 0.8481
AgeGenotypeCloneD 0.0786 0.1643 0.4782 0.6347
AgeGenotypeCloneE 0.0278 0.1643 0.1690 0.8665
AgeGenotypeCloneF -0.0116 0.1643 -0.0705 0.9441
ここで (Intercept) はGenotypeCloneA を基準とした ベースと...
Age はベースとなる傾き
AgeGenotypeClone* は交互作用で，ベースとなる傾きに対する差
ここも参照
また多重回帰分析の場合は，切片以外の項はその変数に対する...
* 文字列ベクトルを数値ベクトルにかえる, 要因をまとめる 2...
as.numeric(as.factor(string.vector))
要因をまとめる crawley, S-Plus, p.301
newgen<- factor(1+(Genotype=="CloneB")+(Genotype=="Clon...
+ 2*(Genotype=="CloneC")+2*(Genotype=="CloneE")
+ 3*(Genotype=="CloneF"))
* lme における anova の解釈 2006 03 17 [#edf6a6ce]
Crawley, S-Plus, p.679, p.690 - 691
anova(k.model1, k.model2 )
# crawley, p.679 複雑なモデル k.model2 が AIC が小さく...
# 説明力( ## Model df AIC BIC logLi...
Test L.Ratio p-value
#k.model1 1 3 2525.512 2537.162 -1259.756 ...
#k.model2 2 4 1932.097 1947.630 -962.048 1 vs 2 59...
しかし AIC が低くても，p-value が大きければ，説明力はない
* [S] S-Plus での画像出力 2006 03 17 [#a03d584c]
http://www.math.aau.dk/~dethlef/Links/latex_figures.html
There are two possibilities for directly producing a
PostScript version of an Splus graphic:
use the menu in the graphics window, with the print
command set to "echo",
so the graph is saved to a file (ps.out.001.ps etc.)
instead of being sent to the printer.
Alternatively, from a program use something like:
postscript(file="foo.ps", horiz=F, onefile=F, print.it=F)
plot...
dev.off()
unix("ps2epsi ./foo.ps ./foo.eps")
The output will be written to the specified file.
It may be useful to process the PostScript file
with the ps2epsi program (supplied with ghostscript)
before including it into latex, since Splus
tends to leave large margins:
ps2epsi foo.ps foo.eps
This also adds a bitmap version of the image,
which is apparently used by certain Macintosh
software to display it on the screen.
* [S] unix への操作 2006 03 17 [#e664fa12]
http://www.uni-muenster.de/ZIV/Mitarbeiter/BennoSueselbec...
# create a character vector from the contents of a file,
# one element of the vector per line of the file
unix("cat myfile")
unix("cat myfile",output=F) # same as
!cat myfile
unix("more", c("line 1", "here is line 2", "and line 3"...
# C shell commands will work if S_SHELL is set to "/bin...
# (You can set S_SHELL explicitly, or let it default to
SHELL or /bin/sh)
!alias # works if S_SHELL is /bin/csh
# You can use unix.shell in place of unix to run csh co...
unix.shell(command = "alias", shell = "/bin/csh", out=F)
# You can make a convenient function to run csh
commands
"csh"<-
function(cmd, ...)
{
unix.shell(command = cmd, shell = "/bin/csh", ....
}
* [S] 余白の大きさを変える 2006 03 17 [#gaa6256e]
以下からの引用です！
http://www.msi.co.jp/splus/support/salon/mcourse/mc12.html
par(mar=c(2, 3, 2, 2))
# 余白の大きさを，行(フォント)の高さを単位に変更
などです。
これらパラメータを設定されると、次に par で同じパラメー...
まで、 あるいは、グラフ・ウィンドウが閉じられるまで継...
忘れないでください。
また S-plus 日本語訳の p.81
ただし lattice パッケージの trellis.device については
ligges, p.167
* [S] Java アプレットでのオンラインヘルプ 2006 03 18 [#...
S-plus のコマンド状で
help.start()
Starting up Java help system. This may take a minute.
終えるには
help.off()
* Regression constant (切片)を除く指定と，その意味 2006 ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.132, p.136
切片(総平均)を除くにはモデル構造式に -1 を含めるが，その...
* 要因を数値に変える unclass 2006 03 18 [#vcc6997d]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* 様々な Anova I, II, III 2006 03 20 [#ra068f7c]
''以下は 次からの引用です！ 石田''
http://www.aichi-gakuin.ac.jp/~chino/anova/chapter1/sec1-...
他に [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.135
SAS では、glm プロシジャで、分析に際して４つの 平方和が用...
その選択には 注意を要する。
まず、Type I 平方和は、逐次平方和 (sequential sums of squ...
(主として釣り合い型デザインの場合）には、この平方和は投入...
Type II 平方和は、当該要因と交互作用を持つような要因や同...
Searle (1987) によれば、Type III 平方和は、Σ 制約付 モデ...
Type III 平方和も Type IV 平方和も、偏平方和と呼ばれるこ...
釣り合い型デザインでは、４つの平方和は、すべて等しい。ま...
* Jittering で重複点をずらしてプロットする 2006 03 20 [#...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.139
points(jitter(Fcat[partner.status == 'low']),
conformity[partner.status == 'low'], pch = 'L')
points(jitter(Fcat[partner.status == 'high']),
conformity[partner.status == 'high'], pch = 'H')
* 要因Aと他要因 BCD全体との差，BCD 内の差を一度に検定 20...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.142 - 143.
* 欠損値の扱い 2006 03 21 [#s6421c0b]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.151
na.action(引数) で設定するが，デフォルトは
R は na.omit (欠損値を含まないデータで解析)
S は na.fail (エラーを出して解析を中止)
S では，計算の際に引数として na.action = na.omit, na.act...
が必要。この二つの違いについては [[John Fox:http://socser...
* offset の意味 2006 03 21 [#h6f113a3]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.153
R のみでのオプション
線形モデルでは，目的変数からその変数を引くのに等しい。
Faraway, GLM, p.63
rate mode で使う
* 対話的にテキストを挿入 locator 2006 03 22 [#if668996]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.165
text(locator(2), c("Close", "One-Sided"))
* contr.poly の出力の L, Q について 2006 03 22 [#ua736596]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.130
* car ライブラリの expand.grid 2006 03 23 [#n73552d4]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.170, p.181
データフレームの行数に合わせて，各列の変数を繰り返す
* [S] テーブルをデータフレームへ変換 2006 03 24 [#r23121...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.184
cbind(expand.grid(dimnames(table)))
Freq = as.vector(table)
* logit と probit の違い 2006 03 28 [#sb64a734]
次の本およびサイトを参照しました．
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
http://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%83%83%...
Logit(p) = log(p/(1-p)) = log(p) - log(1 - p)
''以下は 次からの引用です！ 石田''
http://www.geocities.com/hikachu/2BAStatar.html
Logit&Probit
従属変数が0-1の値をとる場合の回帰分析です．そのまま通常の...
そうすれば予測値は見事[0,1]に収まるというわけです．
この変換にロジスティック関数を用いるのがLogit，累積正規分...
本当はNormitの方が相応しい名前かもしれませんが…．LogitとP...
若干の違いがあるとすれば 1.レアなイベント 2.拡張版の相性 ...
* drop1 は説明変数単独の有意度を測る 2006 03 28 [#ue93c5...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
full モデルとの差を測る
* 信頼区間を求める confint 2006 03 28 [#z9e323cc]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* R with 関数 2006 03 30 [#u4cb5e2d]
以下はRjpwikiより
データを引数 data = hoge と指定できない場合などに便利
例えば interaction.plot を attach されていないデータ trou...
with(troutegg, interaction.plot(period, location, elogi...
を実行する。
interaction.plot(period, location, elogits, data = trou...
* ケースとコントロール 2006 03 31 [#p2df496e]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
終了行:
個人メモを整理中です．不適切な記述が多々あるかと思います...
#contents
* 互いに相関しているサンプルの分散の差の検定 2006 01 06 ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* 偏相関係数の計算方法 2006 01 06 [#w512ce2d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
lm で，対象としたい変数を取り除いたモデルを update で生成...
その変数の 偏 SSR (regression sum of squares) が求まる。
これを SST で割り，平方に開けば良い。
* SST (total sum of squares) の計算方法 2006 01 06 [#z45...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
lm(X ~ 1)
* Anova 表の見方，RSS 残差平方和 2006 01 06 [#xe0ae3d7]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
RSS は回帰によって説明されない変動，長谷川「多変量」 p.22
* [S] S-Plus で locale のエラー [#bb3c3fa8]
Warning: Locale not supported for XmbTextListToTextProp...
Warning: Cannot convert XmString to compound text
http://s.isac.co.jp/FAQ/faq-0032.html を参照
http://S.isac.co.jp/docs-4.5/machdep_3.htmlを参照
その他，myTurboAMD.html の S-plus インストールを参照
* [S] detach rm の順番 2006 01 09 [#qd4b3d6a]
まず detach し，次に rm する。逆に行なうと，
Problem: Object "factories" not found
Use traceback() to see the call stack
とエラー表示
* ネストの式 2006 01 11 [#b93e19cd]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
y ~ A/B
y ~ A * A:B
y ~ A + B %in % A
に等しい。
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.135 - 136
複数の要因からなるデータの nest は，ある要因の水準ごとに
別の要因の傾きを算出する。そして
A) すべての切片が 0 か
B) ある要因のある水準に対する別の要因の傾きがすべて 0 か
を検定する。
ネストモデルについては三中先生のサイト
http://cse.niaes.affrc.go.jp/minaka/R/NestedANOVA.html
を参照
* 切片を除く．その解釈 2006 01 11 [#t4e8647c]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p. 132 下, ...
y ~ x - 1
は説明変数が連続量の場合，切片を 0 にし，全体平均と 全体 ...
すべての説明変数がカテゴリデータの場合，各レベルの平均を...
( -1 が無ければ，全平均，また平均からの差になる)。
* scale location plot 2006 01 14 [#n7eb0bef]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
回帰分析の結果を plot の引数に与えた場合の第三のプロット...
分散が定数であるかを確かめるのに有効。
* カテゴリベクトルの順序づけ ordered 2006 01 14 [#f2988f...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.130
* 変数削減の方法 anova と summary の違い 2006 01 18 [#ad...
Dalgaard, S-Plus, p.181, p.155
summary() で表示された P 値は，表示順に関係なく削除可能。
これに対して anova() の場合は，表示順序が重要。
また高次の相互作用をチェックする場合，ある単独ファクター...
それを含む高次作用が有意の場合は，その単独ファクターは残...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* 適合値と分散の対応 2006 01 18 [#f519a964]
要因A と要因Bの各レベル語とに組合せにおける実験回数が少な...
多少，適合値に対して分散が減少する傾向があっても問題ない。
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* [S] AIC の計算方法 2006 01 27 [#h08b644e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
mode.1<- lm(growth ~ 1)
AIC(model)
* NULL model はパラメータを grand mean 以外推測しない 20...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
model<- aov(y ~ 1)
model<- lm(y ~ 1)
* Helmert contrasts vs. Treatment contrasts 2006 01 30 [...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
y = u + b1 x1 + b2 x2 + b3 x3 + b4 x4
というモデルが想定される場合，五つのパラメータを想定する...
treatment mean は
b1 を 0 と見なし，x1 の平均を u とし，他の b は，u との差...
これに対して Helmert Contrasts は，
全体平均を最初のパラメータとし
最初と二つ目の平均の平均と，最初の平均との差を第二のパラ...
最初と二つ目の平均の平均と，最初と二つ目さらに三つ目の平...
最初と二つ目さらに三つ目の平均の平均と，全体平均の差を第...
またコントラストを手作業で設定するのは [[John Fox:http://...
他にもここを参照
* 繰り返しのない二元配置分散分析 2006 01 31 [#af9d9e29]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Block, Plot, Split-plot の違い 2006 02 01 [#w63cb6bb]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* tapply と list の組み合わせ 2006 02 01 [#fcee3a67]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
attach(splityield)
names(splityiled)
> names(splityield)
[1] "yield" "block" "irrigation" "density"
"fertilizer"
だとして
tapply(yield, list(block, irrigation), sum)
等のように使う。
list については他に crawley, S-Plus, p.406, p.414
* Interaction が有意の時，main effect は削除できない 200...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Fixed Effects はそのファクターレベルに情報がある 2006 ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* [S] interaction.plot の使い方 2006 02 02 [#l37c66ab]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* [S] model.table の使い方 2006 02 02 [#s2e21fb1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Block と Plot, Split-Plot 2006 02 03 [#qc63bace]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
Block (= Plot) をさらに半分にして Split-Plot と呼ぶ。
* Treatment structure, Error structure 2006 02 03 [#efa1...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* Predict と estimate の違い 2006 02 03 [#u93595c1]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
また fitted と predict の違いは [[Faraway:http://www.amaz...
* Corrected, uncorrected Sum os Squares 2006 02 03 [#g33...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
ある変数の値 に その変数の値 をかけた値の総和 uncorrected...
から
その変数の総和の二乗をその変数の数で割ったものを引くと co...
Corrected, uncorrected Sum os Squares の違いは，前者は変...
* 変数変換について 2006 02 07 [#a1daf86a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* predict の使い方 2006 02 06 [#ra24fe2d]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* sample with replacement 2006 02 08 [#v2aabfa4]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
sample with replacement, 取り出した石をもとに戻す
sample without replacement もとに戻さない Hypergeometric
* ヒストグラムで整数値をバーの中心にする 2006 02 09 [#d6...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
break points を 0.5 間隔にする。
hist(mass, breaks = -0.5:16.5)
* Deviance (乖離度) について 2006 02 11 [#x276d817]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
単独モデルの場合
total deviance と residual deviance は
それぞれ satured model との差を表している。
複数モデル比較の場合
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
「パラメータの少ない方のモデルがフィットしている」が帰...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
regression 等ではモデルの当てはめ具合を SSE で測ったが，
GLM では，binomial な proportion データについて deviance ...
また p.526
Deviance is equivalent to sums of squares in linear mod...
p.539 には
total deviance と residual deviance の違い。
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
に単独のモデルについての deviance と，
二つの入れ子になったモデル(どちらかがパラメータが多い)を...
* GLM 小サンプルの weight を小さくする 2006 02 14 [#n8ef...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
fail <- total - success
y <- cbind(success, fail)
model1 <- glm(y ~ explanatory,binomial)
* カテゴリ変数のレベル変更 (アルファベット順以外) 2006 0...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
ordered(target, levels = c("E","D","C","B","A"))
* binomial, exponential, poisson エラーについて(overdisp...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
Everit & Hothorn p.103 - 105.
dispersion parameter は 1 にセットされる。
なぜなら binomial, exponential, poisson では分散は，デー...
residual deviance が residual degrees of freedom よりも大...
overdispersion
この時，モデルの比較はデビアンスではなく，F検定量で行う (...
dispersion パラメータの推定方法は Faraway, p.47, p.60
* logistics analysis of covariance with proportion data ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
例えば三つの薬品のそれぞれについて，分量を変えながら，効...
* wilcox.testが有効ではないケース 2006 02 17 [#q7dce0b3]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
0,1 からのみなるようなデータの場合同順位が多くなる。
* G test，chisq とは別の方法による独立性の検定 2006 02 1...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
glm.crss<- glm(c(22,61,69) ~1, family=poisson)
1- pchisq( glm.crss[13]$deviance, glm.crss[7]$df.residu...
あるいは
glm.crss<- glm( c(15,10,10,15) ~ 1, poisson)
summary( glm.crss)
1 - pchisq(glm.crss$deviance, glm.crss$df.residual)
[1] 0.5695988]
これは
pchisq(glm.crss$deviance, glm.crss$df.residual,
lower = FALSE)
と同じこと
* nuisance 変数 2006 02 18 [#gb4dd871]
crawley, p.551
* 高次の交互作用がある場合，それより低次の要因は省けない...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
* モデル比較を多数くり返す場合の p 値 2006 02 18 [#n9a36...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
"With this many comparisons, I would always work at
p = 0.01 for significance, rather than p = 0.05 "
* 二つの要因に有意差があると出力されても，それは Interce...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
summary(p.580.model)
Call: glm(formula = cbind(relief, patients - relief)
~ treat, family = binomial)
Coefficients:
Value Std. Error t value
(Intercept) -1.067841 0.2471428 -4.320744
treatB 1.959839 0.3427433 5.718094
treatC 2.468734 0.3666004 6.734128
# これは Intercept と差
* GAM (Generalised Additive Models 2006 02 24 [#x27432da]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
y ~ s(x) + s(y) + s(z)
y ~ s + lo(w) + z
s は s smoother , lo は lo smoother を表す。
* Trend がある 2006 02 28 [#v681bb1e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
回帰分析をして，傾きに有意差があれば，そこにトレンドが認...
* プロット描画領域にテキストを加える 2006 03 02 [#wa2f9c...
text 関数を使う。例えば
text(grexp.obs$N[nrow(grexp.obs) ] ,
grexp.obs$K[nrow(grexp.obs) ], labels = "grexp")
これはデータフレーム grexp.obs の N 列の最大値を X 軸に，...
baayen.R を参照。
* スケールの違うプロットを重ねる 2006 03 03 [#p0db652a]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
par(mar = c(5,4,4,5) + 0.1) としてマージンを広げ
tsplot(cinnabar, ylab = "cinnabar") 最初のグラフを描き
par(new = T, xaxs = "d") として新しいプロットが先のプロッ...
* 空のリストを作成する方法 2006 03 03 [#n3bb2f43]
Ligges p. 103
X <- vector(mode = "list", length=10)
* リストに対する演算 2006 03 03 [#w35e23a7]
X <- vector(mode = "list", length=5)
# 5個の要素からなるリストを作成し
lapply(lapply(compare.K, "[", , 2 )
lapply(lapply(compare.K, "[", , "N" )
# リストの要素(データフレーム)のすべてから，
# すべての行の 2 列目(あるいはラベル名が N の列)を取り...
lapply(lapply(compare.K, "[", , 2 ) , "max")
# リストから上記の条件で抽出した各要素について，
# その最大値を取り出す
max(as.numeric (lapply(lapply(compare.K, "[", , 1 ) ,
"max") ))
# 最終的に，リスト全要素の 2 列目の最大値を取り出す
* sapply の使い方 2006 03 04 [#k8e3c932]
Z.list の各要素はベクトルだとする。
unlist(Z.list) でベクトルが返される。
ちなみに sapply(Z.list, "[",1)
は最初の要素をベクトルとして取り出す。
sapply(Z.list, "[") はリストの各要素(ベクトル)を列とした...
Z.data<- data.frame(sapply(Z.list, "["))
とすればデータフレームに変更可能
* paste, parse, expression, eval によるデータフレームの...
x <- "col"
y <- 1:5
z <- "th = 0"
xyz <- paste(x,y,z,sep="")
# xyz
# [1] "col1th = 0" "col2th = 0" "col3th = 0"
"col4th = 0" "col5th = 0"
xyz <- paste(xyz, collapse=",")
# xyz
# [1] "col1th = 0,col2th = 0,col3th = 0,
col4th = 0,col5th = 0"
d <- "data.frame("
d1 <- ")"
d2 <- paste(d,xyz,d1)
# d2
# [1] "data.frame( col1th = 0,col2th = 0,col3th = 0,
col4th = 0,col5th = 0 )"
parse(text =d2)
# parse(text =d2)
# expression(data.frame(col1th = 0, col2th = 0,
col3th = 0, col4th = 0,
# col5th = 0))
Z <- eval(parse(text =d2) )
# Z <- eval(parse(text =d2) )
# Z
# col1th col2th col3th col4th col5th
# 1 0 0 0 0 0
* 文字列の処理，文字列を式に変更 2006 03 04 [#s9534d8b]
Ligges, p.54, p.91
例えばデータフレームの列名をベクトルに変換するには
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
pig<- read.table(
"/S/crawley/data/pig.txt",header=T)
attach( pig )
names( pig )
deparse( names( pig ))
# "c(\"Pig\", \"t1\", \"t2\", \"t3\", \"t4\", \"t5\",
\"t6\", \"t7\", \"t8\", \"t9\")"
parse(text = deparse( names( pig )))
#expression(c("Pig", "t1", "t2", "t3", "t4", "t5",
"t6", "t7", "t8", "t9"))
eval(parse(text = deparse( names( pig ))))
#[1] "Pig" "t1" "t2" "t3" "t4" "t5"
"t6" "t7" "t8" "t9"
parse(text = eval(parse(text = deparse( names( pig )))))
#expression(Pig, t1, t2, t3, t4, t5, t6, t7, t8, t9)
よって次のように行なう
pig.colnames < - eval(parse(text = deparse(names(pig))))
* デ−タフレ−ムやリストをベクトルに変換 2006 03 06 [#n4d2...
x1<- data.frame(y1 = c(1,2,3), y2 = c(11,22,33))
x2<- data.frame(y1 = c(4,5,6), y2 = c(44,55,66))
x3<- data.frame(y1 = c(7,8,9), y2 = c(77,88,99))
X123<- list(x1,x2,x3)
リスト X123 の各要素はベクトルだとする。
unlist(X123) でベクトルが返される。
# データフレームの各列をつなげる
unlist(x1)
あるいは
attach(x1)
z<- NULL
for(i in 1:2){
z<- c(z, get(paste("y", i ,sep="")))
}
# 列数が少ないなら
c(y1,y2) #でも OK
* rep()関数の使い方 2006 03 06 [#w4f9f6e7]
> rep(c("Carroll","Dickens", "Eliot", "Bronte",
"Stevenson", "Gaskell"), rep(3,6))
あるいは
> rep(c("Carroll","Dickens", "Eliot", "Bronte",
"Stevenson", "Gaskell"),each = 3)
* 文字列操作 2006 03 06 [#m580f36e]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-boo...
Ligges p.54
例えば strtrim(labs,2) は文字列ベクトル labs のそれぞれの...
* アルファベットの略語 2006 03 07 [#w7fcf51b]
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
abbreviate(str)
* 対数尤度比を算出 ，松原「統計学の考え方」 2006 03 10 [...
松原 「統計学の考え方」, p82 - p.84
> ftable(p.73)
death Y N
suspect victim
Wh wH 19 132
bL 0 9
Bl wH 11 52
bL 6 97
m2 iterations: deviation 8.881784e-16
$lrt
[1] 0
$pearson
[1] 0
$df
[1] 0
$margin
$margin[[1]]
[1] 1
$margin[[2]]
[1] 2
$margin[[3]]
[1] 1 2
loglin(matu, c(1,2), param =T)
2 iterations: deviation 0
$lrt
[1] 0.1658080
$pearson
[1] 0.1690821
$df
[1] 1
あるいは
library(MASS)
matu<- array(data=c(4,6,20,40),
dim=c(2,2), dimnames = list(Row =c("a1","a2"),
Col = c("b1", "b2")) )
loglm(~Row + Col, data = matu)
x<- array(data=c(15, 10, 10, 15), dim=c(2,2))
# = x<- matrix(c(4,6,20,40), nrow = 2)
kekka<- loglin(x, c(1,2), param =T)
1 -pchisq(kekka$lrt,kekka$df)# 対数尤度比 G による検定...
1 -pchisq(kekka$pearson,kekka$df)# カイ二乗検定による検...
* R での時系列分析 2006 03 10 [#f9c56cde]
R-wiki
R の基本パッケージ base, stats 中の時系列オブジェクトの...
時系列データに関する最も基本的な統計量は
自己共分散(auto-covariance)と
自己相関係数(auto-correlation)である。
R の関連する関数は acf(), pacf(), ccf() である。
関数 acf() は時系列オブジェクトの自己共分散と自己相関係...
既定でそれをプロットする。
関数 pacf() は 偏自己相関係数(partial autocorrelations)...
関数 ccf() は二つの一次元時系列間のクロス相関係数
(cross-correlation)と
クロス共分散(cross-cavarinace)を計算する。
acf(x, lag.max = NULL,
type = c("correlation", "covariance", "partial"),
plot = TRUE, na.action = na.fail,
demean = TRUE, ...)
pacf(x, lag.max = NULL, plot = TRUE,
na.action = na.fail, ...)
ccf(x, y, lag.max = NULL,
type = c("correlation", "covariance"),
plot = TRUE, na.action = na.fail, ...)
* 回帰分析の summary の解釈 2006 03 14 [#c2d6ebc2]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.134
summary(prestige.mod.3)
Call:
lm(formula = prestige ~ income + education + type
+ income:type + education:type)
Residuals:
Min 1Q Median 3Q Max
-13.462 -4.225 1.346 3.826 19.631
Coefficients:
Estimate Std. Error t value Pr(>|t...
(Intercept) 2.276e+00 7.057e+00 0.323 0.74...
income 3.522e-03 5.563e-04 6.332 9.62e-...
education 1.713e+00 9.572e-01 1.790 0.07...
typeprof 1.535e+01 1.372e+01 1.119 0.26...
typewc -3.354e+01 1.765e+01 -1.900 0.06...
income:typeprof -2.903e-03 5.989e-04 -4.847 5.28e-...
income:typewc -2.072e-03 8.940e-04 -2.318 0.02...
education:typeprof 1.388e+00 1.289e+00 1.077 0.28...
education:typewc 4.291e+00 1.757e+00 2.442 0.01...
例えば，ここでは，Intercept は type要因の bc 水準で，これ...
income, education は type要因の bc 水準の時の傾き
typeprf, typewc は切片との差
income:typeprof は，type要因の prof 水準の時の傾き
Dalgaard, p.178
Call:
lm(formula = log10(diameter) ~ log10(conc) * glucose,
data = hellung)
Residuals:
Min 1Q Median 3Q Max
-2.672e-02 -4.888e-03 5.598e-05 3.767e-03 1.761e-02
Coefficients:
Estimate Std. Error t value Pr(>|t...
(Intercept) 1.627926 0.033754 48.230 < 2e-...
log10(conc) -0.046716 0.006846 -6.823 1.51e-...
glucose 0.003418 0.023695 0.144 0.8...
log10(conc):glucose -0.006480 0.004821 -1.344 0.1...
まず最初の二つは glucose があるデータの切片と傾き，次の二...
二つの切片の差は 0.003418 , 傾きの差は -0.006480
差は二つとも有意でない(0.88, 0.185)
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
options(contrasts=c("contr.treatment","contr.poly"))
Call: lm(formula = Weight ~ Age * Genotype)
Residuals:
Min 1Q Median 3Q Max
-0.7787 -0.3986 -0.01139 0.4055 0.8445
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 7.5407 0.3853 19.5724 0.0000
Age 0.2931 0.1162 2.5230 0.0150
GenotypeCloneB 1.0401 0.5449 1.9090 0.0623
GenotypeCloneC -0.9487 0.5449 -1.7411 0.0881
GenotypeCloneD 0.5883 0.5449 1.0797 0.2857
GenotypeCloneE -0.9587 0.5449 -1.7596 0.0849
GenotypeCloneF 1.5693 0.5449 2.8802 0.0059
AgeGenotypeCloneB -0.0241 0.1643 -0.1468 0.8839
AgeGenotypeCloneC -0.0316 0.1643 -0.1926 0.8481
AgeGenotypeCloneD 0.0786 0.1643 0.4782 0.6347
AgeGenotypeCloneE 0.0278 0.1643 0.1690 0.8665
AgeGenotypeCloneF -0.0116 0.1643 -0.0705 0.9441
ここで (Intercept) はGenotypeCloneA を基準とした ベースと...
Age はベースとなる傾き
AgeGenotypeClone* は交互作用で，ベースとなる傾きに対する差
ここも参照
また多重回帰分析の場合は，切片以外の項はその変数に対する...
* 文字列ベクトルを数値ベクトルにかえる, 要因をまとめる 2...
as.numeric(as.factor(string.vector))
要因をまとめる crawley, S-Plus, p.301
newgen<- factor(1+(Genotype=="CloneB")+(Genotype=="Clon...
+ 2*(Genotype=="CloneC")+2*(Genotype=="CloneE")
+ 3*(Genotype=="CloneF"))
* lme における anova の解釈 2006 03 17 [#edf6a6ce]
Crawley, S-Plus, p.679, p.690 - 691
anova(k.model1, k.model2 )
# crawley, p.679 複雑なモデル k.model2 が AIC が小さく...
# 説明力( ## Model df AIC BIC logLi...
Test L.Ratio p-value
#k.model1 1 3 2525.512 2537.162 -1259.756 ...
#k.model2 2 4 1932.097 1947.630 -962.048 1 vs 2 59...
しかし AIC が低くても，p-value が大きければ，説明力はない
* [S] S-Plus での画像出力 2006 03 17 [#a03d584c]
http://www.math.aau.dk/~dethlef/Links/latex_figures.html
There are two possibilities for directly producing a
PostScript version of an Splus graphic:
use the menu in the graphics window, with the print
command set to "echo",
so the graph is saved to a file (ps.out.001.ps etc.)
instead of being sent to the printer.
Alternatively, from a program use something like:
postscript(file="foo.ps", horiz=F, onefile=F, print.it=F)
plot...
dev.off()
unix("ps2epsi ./foo.ps ./foo.eps")
The output will be written to the specified file.
It may be useful to process the PostScript file
with the ps2epsi program (supplied with ghostscript)
before including it into latex, since Splus
tends to leave large margins:
ps2epsi foo.ps foo.eps
This also adds a bitmap version of the image,
which is apparently used by certain Macintosh
software to display it on the screen.
* [S] unix への操作 2006 03 17 [#e664fa12]
http://www.uni-muenster.de/ZIV/Mitarbeiter/BennoSueselbec...
# create a character vector from the contents of a file,
# one element of the vector per line of the file
unix("cat myfile")
unix("cat myfile",output=F) # same as
!cat myfile
unix("more", c("line 1", "here is line 2", "and line 3"...
# C shell commands will work if S_SHELL is set to "/bin...
# (You can set S_SHELL explicitly, or let it default to
SHELL or /bin/sh)
!alias # works if S_SHELL is /bin/csh
# You can use unix.shell in place of unix to run csh co...
unix.shell(command = "alias", shell = "/bin/csh", out=F)
# You can make a convenient function to run csh
commands
"csh"<-
function(cmd, ...)
{
unix.shell(command = cmd, shell = "/bin/csh", ....
}
* [S] 余白の大きさを変える 2006 03 17 [#gaa6256e]
以下からの引用です！
http://www.msi.co.jp/splus/support/salon/mcourse/mc12.html
par(mar=c(2, 3, 2, 2))
# 余白の大きさを，行(フォント)の高さを単位に変更
などです。
これらパラメータを設定されると、次に par で同じパラメー...
まで、 あるいは、グラフ・ウィンドウが閉じられるまで継...
忘れないでください。
また S-plus 日本語訳の p.81
ただし lattice パッケージの trellis.device については
ligges, p.167
* [S] Java アプレットでのオンラインヘルプ 2006 03 18 [#...
S-plus のコマンド状で
help.start()
Starting up Java help system. This may take a minute.
終えるには
help.off()
* Regression constant (切片)を除く指定と，その意味 2006 ...
[[Crawley S-PLUS:http://www.bio.ic.ac.uk/research/mjcraw/...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.132, p.136
切片(総平均)を除くにはモデル構造式に -1 を含めるが，その...
* 要因を数値に変える unclass 2006 03 18 [#vcc6997d]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* 様々な Anova I, II, III 2006 03 20 [#ra068f7c]
''以下は 次からの引用です！ 石田''
http://www.aichi-gakuin.ac.jp/~chino/anova/chapter1/sec1-...
他に [[John Fox:http://socserv.mcmaster.ca/jfox/]] p.135
SAS では、glm プロシジャで、分析に際して４つの 平方和が用...
その選択には 注意を要する。
まず、Type I 平方和は、逐次平方和 (sequential sums of squ...
(主として釣り合い型デザインの場合）には、この平方和は投入...
Type II 平方和は、当該要因と交互作用を持つような要因や同...
Searle (1987) によれば、Type III 平方和は、Σ 制約付 モデ...
Type III 平方和も Type IV 平方和も、偏平方和と呼ばれるこ...
釣り合い型デザインでは、４つの平方和は、すべて等しい。ま...
* Jittering で重複点をずらしてプロットする 2006 03 20 [#...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.139
points(jitter(Fcat[partner.status == 'low']),
conformity[partner.status == 'low'], pch = 'L')
points(jitter(Fcat[partner.status == 'high']),
conformity[partner.status == 'high'], pch = 'H')
* 要因Aと他要因 BCD全体との差，BCD 内の差を一度に検定 20...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.142 - 143.
* 欠損値の扱い 2006 03 21 [#s6421c0b]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.151
na.action(引数) で設定するが，デフォルトは
R は na.omit (欠損値を含まないデータで解析)
S は na.fail (エラーを出して解析を中止)
S では，計算の際に引数として na.action = na.omit, na.act...
が必要。この二つの違いについては [[John Fox:http://socser...
* offset の意味 2006 03 21 [#h6f113a3]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.153
R のみでのオプション
線形モデルでは，目的変数からその変数を引くのに等しい。
Faraway, GLM, p.63
rate mode で使う
* 対話的にテキストを挿入 locator 2006 03 22 [#if668996]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.165
text(locator(2), c("Close", "One-Sided"))
* contr.poly の出力の L, Q について 2006 03 22 [#ua736596]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.130
* car ライブラリの expand.grid 2006 03 23 [#n73552d4]
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.170, p.181
データフレームの行数に合わせて，各列の変数を繰り返す
* [S] テーブルをデータフレームへ変換 2006 03 24 [#r23121...
[[John Fox:http://socserv.mcmaster.ca/jfox/]] p.184
cbind(expand.grid(dimnames(table)))
Freq = as.vector(table)
* logit と probit の違い 2006 03 28 [#sb64a734]
次の本およびサイトを参照しました．
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
http://ja.wikipedia.org/wiki/%E3%83%AD%E3%82%B8%E3%83%83%...
Logit(p) = log(p/(1-p)) = log(p) - log(1 - p)
''以下は 次からの引用です！ 石田''
http://www.geocities.com/hikachu/2BAStatar.html
Logit&Probit
従属変数が0-1の値をとる場合の回帰分析です．そのまま通常の...
そうすれば予測値は見事[0,1]に収まるというわけです．
この変換にロジスティック関数を用いるのがLogit，累積正規分...
本当はNormitの方が相応しい名前かもしれませんが…．LogitとP...
若干の違いがあるとすれば 1.レアなイベント 2.拡張版の相性 ...
* drop1 は説明変数単独の有意度を測る 2006 03 28 [#ue93c5...
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
full モデルとの差を測る
* 信頼区間を求める confint 2006 03 28 [#z9e323cc]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
* R with 関数 2006 03 30 [#u4cb5e2d]
以下はRjpwikiより
データを引数 data = hoge と指定できない場合などに便利
例えば interaction.plot を attach されていないデータ trou...
with(troutegg, interaction.plot(period, location, elogi...
を実行する。
interaction.plot(period, location, elogits, data = trou...
* ケースとコントロール 2006 03 31 [#p2df496e]
[[Faraway:http://www.amazon.co.jp/Extending-Linear-Model-...
ページ名:
