ぎょーむ日誌 2006-04-07
2006 年 04 月 07 日 (金)
-
0820 起床．
朝飯．
コーヒー．
さっさと大学に
……
と出ようとしたところ，
昨日「久保さん，
セーターの肘に穴があいてるじゃないですかっ」
と大学院生に叱られてしまったことを思いだした．
それほどでかい穴ではないので裏側からぬいあわせて処置．
0955 自宅発．
晴．
1010 研究室着．
-
昨日，
粕谷さんからいただいた質問．
{応答した, 応答しなかった} 列のくみあわせを，
{0, 1} vector にして
data.frame()
つくりなおすという問題．
たとえば
c(2, 1) を c(1, 1, 0)
とするような
……
> test1g
x noresp resp total ida
1 1 4 1 5 1
2 2 3 2 5 2
3 3 4 1 5 3
4 3 1 4 5 4
5 4 2 3 5 5
6 5 1 4 5 6
というデータフレームがもとです。説明変数が1のときに反応個体数が1で非反応個体
数が4だった、という内容です。ここからtest1gyという以下のものを作ります。
> test1gy
x noresp resp total ida
1 1 0 1 5 1
2 1 1 0 5 1
3 1 1 0 5 1
4 1 1 0 5 1
5 1 1 0 5 1
6 2 0 1 5 2
7 2 0 1 5 2
8 2 1 0 5 2
9 2 1 0 5 2
10 2 1 0 5 2
... (以下略)
なぜかかる変換が必要かといえば，
我らが親愛なる (そしてややまぬけな)
glmmML()
はこういった {0, 1} 応答変数でなければ
random effects ありの
二項分布モデル (family = binomial)
のパラメーター推定計算やってくれないためだ．
他の GLMM 計算関数
だと cbind(resp, noresp) ~ ...
と書けるのだが．
この質問メイルにつけられていた粕谷さんのコード例は
for (...) { ... }
をたくさん使ったものだったので，
そうではない別の例を作ってみた．
私の返信
(以下で使ってるファイル;
conv01.R,
test1g.csv)．
久保です．R ふうコードというべきかどうかわかりませんが，極端な一例とし
て for を使わない変換関数を作ってみました．
# ------------------------------------------------------------------
conv01 <- function(data, resp, noresp, col.join = NULL)
{
# check errors
if (length(resp) != length(noresp))
stop("ERROR: length(resp) != length(noresp)")
if (length(col.join) < 2 & is.null(col.join))
col.join <- colnames(data) # join all columns
# {resp, noresp} --> {0, 1} vector
resp01 <- c(apply(
data.frame(resp, noresp), 1,
function(x) c(rep(1, x[1]), rep(0, x[2]))
))
# join data columns
total <- resp + noresp
cbind(
data.frame("resp01" = resp01), # {0, 1} vector
data[
sapply(1:nrow(data), function(x) rep(x, total[x])),
col.join
]
)
}
# ------------------------------------------------------------------
これを conv01.R といった名前のテキストファイルに保存しておき，
source("conv01.R") とすると conv01.R をよみこんで関数 conv01() の定義
を記憶します．
データが test1g.csv
# ------------------------------------------------------------------
x,noresp,resp,total,ida
1,4,1,5,1
2,3,2,5,2
3,4,1,5,3
3,1,4,5,4
4,2,3,5,5
5,1,4,5,6
# ------------------------------------------------------------------
だとすると
d <- read.csv("test1g.csv")
として
conv01(d, resp = d$resp, noresp = d$noresp)
だの
conv01(d, resp = d$resp, noresp = d$noresp, col.join = c("x", "ida"))
とすれば変換した data.frame を返してきます．
他にもメイルやりとりとか．
粕谷さんからおススめ本情報
(といっても出版予定は 10 月ってことだが)
Applied Regression and Multilevel/Hierarchical Models (Analytical Methods
for Social Research) (Hardcover)
by Andrew Gelman, Jennifer Hill.
* Publisher: Cambridge University Press (October 31, 2006)
* Language: English
* ISBN: 0521867061
著者の
本ペイジ
に目次などあり．
昼飯．
苫小牧樹木直径成長解析の論文原稿の世界，
にはいりこもうとする努力．
しかしプリンター雑用がわりこんでくる．
しかし新規発注も修理もススまぬ状況，
とわかった．
北大生協の説明．
また、伝票処理についてですが、
調達課から「納品の日付と伝票の日付を同じにする」との通達をいただいております。
発注の手配等をかけることは可能ですが、実際に納品・修理を行うのは、予算が執行されてから、
ということになります。
お待ちいただく形になってしまいますが、ご了承いただきますようよろしくお願いします。
予算の執行っていつになるのやら
……
そして本日は R
こんさるメイルが多い．
夜になってもつづく
……
この季節になんでこんなにデータ解析の問い合わせが増えるのだろうか?
上の「予算の執行」問題，
北大生協からいつになったら金が出せるのかといった質問メイルがきたので，
そのあたりの件について会計係にといあわせのメイルだすと
……
なぜか電話で返事いただき
「本日から公費 (校費?) なら使える」
とのこと．
科研費はとうぜんながらまだ．
さて公費 (校費?)
でしはらってよいのか
……
このあたりの基本がわかってないので，
来週にでも雪野さんにお尋ねするのがよろしかろう．
といったどたばたにまきこまれてるうちに空腹になってきたので撤退．
1845 研究室発．
1900 帰宅．
晩飯．
[今日の運動]
[今日の食卓]
- 朝 (0900):
米麦 0.6 合．
ネギ・しらたき・豆腐の煮物．
- 昼 (1310):
研究室お茶部屋．
米麦 0.6 合．
ネギ・しらたき・豆腐の煮物．
- 晩 (2000):
米麦 0.7 合．
ネギ・しらたき・豆腐の煮物．