Rを用いた「繰り返しがある直交表実験計画法」の分析方法

菅民郎著「Excelで学ぶ実験計画法」を読んで、実験計画法を勉強しています。
本書はExcelで分析しておりますが、Rでも分析できるようになりたいと考えています。

本書7.2節の「繰り返しがある場合の直交表実験計画」で躓きました。内容は、
　・２水準の因子が４つ（因子A, B, C, D）ある実験を想定
　・直交表L8を利用して、組合せ数を16→8に削減
　・被験者3名から各組合せに対する評価を取得
　　　⇒　すなわち、8×3=24データがある
というものです（添付画像が組合せと被験者の評価になります）。

被験者という因子による影響を考慮するため、①組合せ（因子N）と被験者（因子R）を繰り返しのない二元配置分散分析を実施し、②因子A, B, C, Dの分散分析を反映するという流れで分析されていました。

本書で解説されているような、被験者複数名から各組合せに対する評価を取得するような事例が、ネット上を探してみてもなかなか見つかりません。
特にRによる実験計画法に関しては、そもそも情報がほとんどなく（下記スライドのようにデータが１つのものは見つかりますが）、困っております。
　https://www.slideshare.net/itoyan110/r-14261638

被験者複数名から各組合せに対する評価を取得した場合、Rを用いてどのように分析すればよいか、ご教示いただきたいです。
よろしくお願いいたします。

No.1 ベストアンサー 回答者： kamiyasiro 回答日時： 2023/08/01 22:25

L8直交表を生成するステップや、因子を割り付けるステップはどうでもよく、ANOVAを作るだけで良いですか？

具体的には、掲載された表は、1,2,4,7列に割り付けてありますが、ここは既存の表と線点図などを使って手動で行う、すなわち統計パッケージのように自動的に生成されなくても良いのですね？

解析ステップは、
①「被験者間に差があるとは言えない」ことを調べるため二元配置の分散分析を行う
②「被験者は単なる繰り返し」と考えてL8直交表の分散分析を行う

これをRで行いたい、ということですね？

方法は、たぶん、（自分でも暇なときに確認したいと思いますが）

①は、行をｒ、列をｃとすると、
r1,c1,3
r1,c2,4
r1,c3,4
r2,c1,4
r2,c2,3
r2,c3,2
r3,c1,5
r3,c2,5
・
・
・

②は、因子をa,b,c,dとすると、
a1,b1,c1,d1,3
a1,b1,c1,d1,4
a1,b1,c1,d1,4
a1,b1,c2,d2,4
a1,b1,c2,d2,3
a1,b1,c2,d2,2
・
・
・

というデータを作り、最後の列をoutcomeとかいう列名にし、
anova(lm(outcome ~ ., data = x))

とやれば、分散分析表が出てくるハズです。

因子を数値にすると、違った答えが出て来てしまいます。factorにしてやることがポイントです。

数値でやるなら、ー１と１というコード化された状態でないとダメだと思います。今のままでは、各列の内積（積和）は０にはなっていませんから。
（確認した訳ではありませんが・・・）

この回答へのお礼

詳しい解説をありがとうございます。Rで実際に分析してみました。

A, B, C, D, R, X
a1, b1, c1, d1, r1, 3
a1, b1, c2, d2, r1, 4
a1, b2, c1, d2, r1, 5
a1, b2, c2, d1, r1, 5
a2, b1, c1, d2, r1, 8
:

というdatデータ（Rは被験者番号）を用意して、

anova(lm(X~R+A+B+C+D+A:B+A:C,data=dat))

とすると、書籍に記載されている分散分析表に近い結果が得られました。

書籍と異なる点は、①不適合と純誤差は得られない、②寄与率は得られない点でした。
②は、書籍に記載されている式にしたがい、算出することはできますが...。

お礼日時：2023/08/02 12:24

No.8 回答者： kamiyasiro 回答日時： 2023/08/03 11:07

#7です。

＞ 応答曲面法について学習するのとあわせて・・・

応答曲面法のように、より近代的になってくると、実はRの方が良いです。

今は、大量のデータが入手可能なので、２次関数ではなく非線形関数で近似してグネグネの応答を考えるのが主流になって来ました（多目的設計探査で検索してみて下さい）。

・gam()やTps()などの関数で近似できる
・optim()関数で最適点探索が容易にできる
・ライブラリGAなど遺伝的アルゴリズムによる探索も可能

これらが無料でできるって、R恐るべしです。

この回答へのお礼

試しに、Rによる応答局面法について調査してみたところ、情報がかなりあったため、もしやと思っておりました。
有益な情報を沢山、ご教示いただき、本当にありがとうございました！

お礼日時：2023/08/03 18:19

No.7 回答者： kamiyasiro 回答日時： 2023/08/02 16:44

#5です。

大切なことを忘れていました。
自動的なモデル選択は、step(result)で出来るのですが、通常の回帰はこれで良いのです。でも、分散分析表のモデル選択はちょっと違います。

分散分析には、「ヘレディティ heredity」という暗黙のルールがあります。
そのルールとは、

①ストロング・ヘレディティ・・AB交互作用があるときは、モデル中にA主効果もB主効果も残さなければならない。

②ウィーク・ヘレディティ・・AB交互作用があるときは、モデル中にA主効果かB主効果かのどちらかを残さなければならない。

世間的には、①が主流です。
このルールは、お手持ちの教科書に書いてないだろうし、ネットで探しても出てきませんが、StatWorksやJMPなどの統計パッケージは、必ず守っています。

例えば、FinFout法で、A主効果を誤差にプーリングしようとすると、
『要因Aに関わる全ての交互作用も誤差にプーリングしますが、よろしいですか』
と聞いています。

この件や、３水準系の分散分析などのこともあって、私はStatWorksなどの統計パッケージを使うようにしています。

ご注意をお願いします。

例えば、

https://www.jstor.org/stable/41806068

自動翻訳すると「モデルの遺伝タイプの制約」と出てくる箇所です。

この回答へのお礼

このような側面を考慮すると、StatWorksやJMPなどのソフトウェアの使用も検討する必要がありそうですね...。
フリーソフトであること、自分で手法の指定ができることから、長らくRを使用していますが、StatWorksもありますし、応答曲面法について学習するのとあわせて、挑戦してみようと思います。

お礼日時：2023/08/02 17:00

No.6 回答者： kamiyasiro 回答日時： 2023/08/02 15:40

#1・#3です。

純誤差の件、Residuals の Mean Sq の平方根で求まりますが、aov関数なら以下のように、

result <- lm(X~R+A+B+C+D+A:B+A:C,data=dat)
aov(result)

とやれば、Residual standard error が表示されます。

次に寄与率の件、
私が知っている限りでは、パッと求めてくれるRの関数は無いと思います。
Rで式を書いていると面倒だし、人的ミスが起きるリスクもありますから、関数があるとありがたいですね。

３水準系の直交表の解析についてですが、
私がRで直交表の解析をやっていたのは何年も前で、その時にトラップを経験しました。
現在は、新たなライブラリが出ているかもしれません。

でも、その後は「応答曲面法」の方が主流になりました。応答曲面関数は高次交互作用はモデルに取り込まない（＝式の最高次数は２次までで交互作用の自由度は１）のです。私の勤務先はメーカーですが、開発現場では既に３水準系のL27なんて使わないのですよ。割り付けできる列数が少なくて非効率だからです。交互作用が自由度４も奪っていくからです。

それに、応答曲面法は、多目的最適化も可能です。
ご質問者様には、３水準系は軽く読み飛ばして応答曲面法に進まれることをお勧めします。

すると、「D最適計画」とか、新しい計画表に出会うことになると思います。

No.5 回答者： kamiyasiro 回答日時： 2023/08/02 15:00

#2です。

まず、プーリングから。

FinFout法によるモデル選択ではなく、最尤法によるモデル選択になりますが、結果は大きく変わりません。

もし、ご質問者様が#1のコメントのような下記のモデル式を用いられたなら、

result <- lm(X~R+A+B+C+D+A:B+A:C,data=dat)

のように、モデルフィットした結果をresult変数に放り込んでおき、

step(result)

でモデル選択ができます。

No.4 回答者： kamiyasiro 回答日時： 2023/08/01 23:29

#3です。

２水準系の交互作用の自由度は１ですが、３水準系の交互作用の自由度が４というのは、LL、LQ、QL、QQの４項を解いているということです。

ちなみに、

LL交互作用は、X字型の要因効果図になります。
LQ交互作用は、山谷逆転型の要因効果図になります。

L18直交表（タグチ一派が使うヤツ）は、主効果にこの高次交互作用が部分交絡しているので、因子を割り付けていない列の要因効果図に幽霊のような効果が出るのです。連中はそんなこと理解していないですが・・・。

てか、ここまで来ると、完全に学習範囲を逸脱していますね。すみません。

No.3 回答者： kamiyasiro 回答日時： 2023/08/01 23:07

#2です。

学習範囲を逸脱してしまいますが、２水準系がRでできたからと言って、決して同じ調子で３水準系をやらないで下さい。

１次モデルをL、２次モデルをQとすると、直交表の解析は、

・２水準系交互作用は、LL項しかない。
・３水準系交互作用は、LL、LQ、QL、QQという高次交互作用項まで含んでいる（十把一絡げで解いている）。

一方、anova(lm(・・・
というRの分散分析表の作り方は、交互作用は例え３水準系であってもLL項までしか解けないのです。

Rでやるときのトラップですので、ご注意を。

この回答へのお礼

トラップに嵌るところでした。細かい部分まで解説いただき、誠にありがとうございます。

書籍でも、３水準系の解説は記載されていますが、交互作用までは言及されていませんでした。

３水準系以上で高次交互作用を取り扱いたい場合は、専門のソフトウェアが必要ということなのでしょうか...。
（今のところ、そういった実験を検討しているわけではありませんが、後学のため）

お礼日時：2023/08/02 12:24

No.2 回答者： kamiyasiro 回答日時： 2023/08/01 22:54

#1です。

①は交互作用が心配されても求められませんが、（繰り返しの無い二元配置計画は交互作用を求めることができないから）
②は交互作用を求めることができます。

・交互作用を求める場合は、２水準系の場合は、データを変更するか、lmの()の中のモデル式を変更すれば可能です。
ANOVA中に交互作用効果も表示することは必要ですか？
（このご質問とは関係ありませんが、３水準系は、簡単な解析方法が思い付きません。なぜならa^2b列とab^2列の２列を跨いで交互作用が出るからです。）

・あと、効果の無い因子を誤差にプーリングしてANOVAを作る必要はありますか？
自動で行う必要はありますか？その場合、RではFinFoutではなく、最尤法になってしまうかも。

この回答へのお礼

②の交互作用も得られました。ありがとうございます。

プーリングについてもできるようになりたいと考えております。
可能であれば、Rによる手法をご教示いただければ幸いです。

お礼日時：2023/08/02 12:24