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

No.17068　RとSASの予測値が不一致　　【赤岳】　2012/06/13(Wed) 13:22

いつもご教示ありがとうございます。
さて，またRで悩みが発生しました。

Rの結果とSASの結果が一部異なり，困っています。
ご教示いただけませんでしょうか。
まず，検討したデータを紹介します。
2種類の薬剤のAUCデータを，9例の被験者を用いて2x2クロスオーバー法で取得しました。
その結果は，次のとおりとなりました。
> Data
Subject Treatment Group Seq AUC
1 N1 C A p2 1.1058813
2 N1 F A p1 1.1149844
3 N2 C A p2 0.8242625
4 N2 F A p1 0.9816918
5 N3 C A p2 0.7548069
6 N3 F A p1 0.9060925
7 N4 C B p1 0.9809392
8 N4 F B p2 0.9104003
9 N5 C B p1 0.8964931
10 N5 F B p2 0.8811163
11 N6 C A p2 0.9328288
12 N6 F A p1 0.7908479
13 N7 C A p2 0.9632446
14 N7 F A p1 0.7843819
15 N8 C B p1 0.8812363
16 N8 F B p2 0.9370914
17 N9 C B p1 0.7809003
18 N9 F B p2 0.8588198

Subjectは被験者番号，Treatmentは薬剤（CとF），Groupは9例を2群に分けたもの，Seqはp1が第1期投与，p2が第2期投与で，p1とp2の間は4週間のWashoutを置いている。
この結果を，混合効果モデルにて，次のモデル式をおいて，Rでプログラミングしました。
Y＝（μ＋δ＋ε）＋β1*Treatment＋β2*Group + β3* Seq
μは平均，δはランダム切片誤差，εは誤差，β1-3はそれぞれの変数の傾き
Subjectを変量効果とおいています。

> library(nlme)
> model.AUC<-lme(AUC ~ Treatment + Group + Seq,random =~1|Subject, method ="REML", data=Data)
> summary(model.AUC)
Linear mixed-effects model fit by REML
Data: Data
AIC BIC logLik
-3.236242 0.598102 7.618121

Random effects:
（省略）

Fixed effects: AUC ~ Treatment + Group + Seq
Value Std.Error DF t-value p-value
(Intercept) 0.9099199 0.05180509 7 17.564296 0.0000
TreatmentF 0.0056798 0.04275350 7 0.132850 0.8981
GroupB -0.0250277 0.06126703 7 -0.408501 0.6951
Seqp2 0.0062849 0.04275350 7 0.147004 0.8873
(省略）

TreatmentFのValueは，0.0056，SE0.00427という結果が得られました。
これは，薬剤CとFの差です。
この結果は，SASの結果とぴったり一致しました。

しかし，Rのデフォルトでは，薬剤CとFそれぞれの予測値が表示されていないので，次のようにプログラミングして，予測値の平均とSDを算出しました。

> Estimate<-predict(model.AUC) #AUCの各例の予測値
> Data.new<-data.frame(Data, Estimate)
> Mean<-aggregate(Data.new$Estimate,list(treatment=Data.new$Estimate, list(treatment= Data.new$Treatment), mean)
> Sd<-aggregate(Data.new$Estimate,list(treatment=Data.new$Estimate, list(treatment= Data.new$Treatment), sd)
> Mean
treatment x
1 C 0.9022881
2 F 0.9072696
> Sd
treatment x
1 C 0.04683475
2 F 0.04493016

しかし，この結果は，SASと一致しませんでした。
SASでは，
Treatment Estimate Se Df
C 0.9005 0.03735 12.5
F 0.9062 0.03735 12.5

このRでの予測値の平均，SdとSASでのそれらとがどうして一致しないのか悩んでいます。
どういう原因が考えられるか，ご教示いただけませんでしょうか。

No.17071　Re: RとSASの予測値が不一致　　【赤岳】　2012/06/14(Thu) 13:01

SASの結果の方が妥当性があるように思います。
Rでの，C群とF群の予測値を引き算すると，0.00498，
SASでの，C群とF群の予測値を引き算すると，0.0057
となり，SASの結果はsummary(model.AUC)のTreatmentFの結果と一致します。
Rの予測値の出し方に問題があるように思っています。
しかし，予測値を出すコマンドは，「predict」か「fitted」しか知りません。
どなたか，Rでの予測値の出し方をご教示いただけませんでしょうか。

No.17072　Re: RとSASの予測値が不一致　　【青木繁伸】　2012/06/14(Thu) 14:37

> TreatmentFのValueは，0.0056，SE 0.00427という結果が得られました。
> これは，薬剤CとFの差です。
> この結果は，SASの結果とぴったり一致しました。

「ぴったり」とはどの程度ですか?
収束計算に使う限界値やアルゴリズムによっても，表示される有効数字以下の違いはあるかもしれません。

> Rでの予測値の出し方をご教示いただけませんでしょうか。

lme が返すオブジェクト（model.AUC）の中に，fitted という要素が既に計算されていますよ（model.AUC$fitted で取り出せる）。

No.17074　Re: RとSASの予測値が不一致　　【赤岳】　2012/06/14(Thu) 16:46

青木先生，
ご連絡ありがとうございます。
SASでもぴったり一致したという点について，
RのSummary：TreatmentF 0.0056798 0.04275350
SASの結果：Treatment 0.00568 0.04275
です。

また，model.AUCfittedで取り出してみました。ご教示ありがとうございます。>Data1<－data.frame(Data,model.AUC$fittedで取り出してみました。ご教示ありがとうございます。 > Data1<-data.frame(Data, model.AUC$fitted)
> Mean<-aggregate(Data1$fixed,list(treatment=Data1$fixed, list(treatment= Data1$Treatment), mean)
> Mean
treatment x
1 C 0.9022881
2 F 0.9072696

この結果は，Rのpredictでの結果と一致しますが，SASとは依然一致しません。
RとSASが結構大きな差があるので，気になっています。

No.17075　Re: RとSASの予測値が不一致　　【青木繁伸】　2012/06/14(Thu) 16:50

> SASでもぴったり一致したという点について，
> RのSummary：TreatmentF 0.0056798 0.04275350
> SASの結果：Treatment 0.00568 0.04275

表示されているよりも下位の桁がどうなっているかがわからないといけません。

No.17076　Re: RとSASの予測値が不一致　　【赤岳】　2012/06/14(Thu) 17:17

入手したSASの結果は下5桁でした。桁数を増やせるようですので，入手してみます（時間がかかるかも知れませんが，ご容赦ください）。