「統計学関連なんでもあり」の過去ログ--- 041

No.06810　lm関数で対比の使い方　　【hiroshi】　2008/06/14(Sat) 22:26

昨日からRを使用し始めたばかりで，どうしても分からないことがあります。
マウスに眼圧負荷(IOP)をかけると網膜神経節細胞(RGC)が減少します。ここに2種類の薬を使用し，網膜神経節細胞の減少に対する効果の違いをみようと思っていますが，実験の前に検出力を確かめたいと思っています。薬1と薬2で各群30匹ずつ使用の予定です。眼圧負荷は0から800まで一様分布します。網膜神経節細胞は無処置眼と処置眼の割合(%)を評価します。薬剤1は RGC = 100 + IOP * -0.07 という式に大体従うことが分かっています。薬剤2は RGC = 100 + IOP * -0.35 という式に従うと予想しています。RGCの測定誤差の標準偏差は約15です。二群で傾きに差があるかどうかを検定したいと思い，以下のようにしました。

library(gmodels)
subjects=30; SD=15; katamuki1=-0.07; katamuki2=katamuki1*0.5;
dat <- data.frame(mouse.id=gl(subjects, 2),
                  kusuri=gl(2, subjects))
dat$IOP <- runif(subjects*2, min=0, max = 800)
dat$RGC <- rnorm(subjects*2, mean=0, sd=SD) + (ifelse(dat$kusuri==1, katamuki1, katamuki2) * dat$IOP) + 100
model <- lm(RGC ~ kusuri:IOP, data=dat)
ans <- fit.contrast(model, "kusuri:IOP", c(1, -1)) #これがうまくいかない

以下のような警告がでてうまくいきません。

Warning message:
In model.matrix.default(mt, mf, contrasts) :
  variable 'kusuri:IOP' is absent, its contrast will be ignored

どのようにすれば宜しいのでしょうか。なにとぞよろしくお願いいたします。

No.06811　Re: lm関数で対比の使い方　　【hiroshi】　2008/06/14(Sat) 22:55

自己レスですみません。
ans <- fit.contrast(model, "kusuri:IOP", c(1, -1))
を以下のようにしたらうまくいった様な気が。
cmat <- c("kusuri1:IOP"=1, "kusuri2:IOP"=-1)
ans <- estimable(model, cmat)
これでいいのでしょうか。全く慣れていないのでとんでもない間違いをしているような気もします。詳しい方よろしくお願いいたします。