中心複合計画と応答曲面法について

統計ソフトにただ数字を入れているだけの統計については全くの素人です。よろしくお願いします。
中心複合計画(Central composite (rotatable) design)を使用して、２個の計画変数（添加濃度Xと反応時間Y）に対する独立変数の応答曲面を出そうとしています。中心4回繰返しで２つのスターポイント(α＝1.4142)をとって、12点の実験を行う計画を立て、結果を出しました。結果はaX^2+bY^2+cXY+dX+eY+fの二次モデルをステップワイズで有意差の無い因子(F<4)を排除しました。

ある結果でステップワイズをかけると完全に式が無くなったり、単純な１次式（dX+fのような）になりましたが、実験している感覚だと三次式的な挙動をしていると感じました。そこで、gX^3+hY^3+(二次モデル）という式を作ってステップワイズをかけたところ高いR^2が出ました。邪道な解析だと思いますが、このような中途半端な三次モデルで現象を説明することは可能でしょうか？
このような方法で論文の考察を行っている例はございますか？

本来は二次モデルで収まるように予備実験すべきなのかもしれませんが、どうしても予備実験も再実験もできない状況でしたので。。。
統計の知識が無いもので、説明自体に誤りがあるかもしれません。何卒よろしくお願いいたします。

No.2 回答者： kamiyasiro 回答日時： 2015/10/19 20:51

#1です。

間違えました。

XTX（Tは転置）を取った時に、右下がり対角線以外が
1になります。
↓
0になります。

No.1 回答者： kamiyasiro 回答日時： 2015/10/18 19:56

企業でSQC推進を担当する統計家です。

まず、どういうソフトを使用して生成されたか知りませんが、
今の計画は、2変数の中心複合計画としては、良くない計画です。
本来、中心点の繰り返しを8回行わないと直交しません。
直交するときは、因子の行列をXとして、
XTX（Tは転置）を取った時に、右下がり対角線以外が
1になります。
現在の計画のXTXは下記のとおりであり、
2次項間に交絡が出ています。
そのために、2次項が打ち消しあっているかもしれません。
（^2は.2と表現）

x1 x2 x1.2 x2.2 x1x2
x1 8 0 0.000 0.000 0
x2 0 8 0.000 0.000 0
x1.2 0 0 6.666 -1.333 0
x2.2 0 0 -1.333 6.666 0
x1x2 0 0 0.000 0.000 4

次に、3次のモデルですが、
中心複合計画のスター点距離を1ではなく、今回のように1.4などの値とすると
軸上に5水準できますので、1見4次のモデル推定まで可能なように見えます。
しかし、高次のモデル
1st・・・x1，x2
2nd・・・x1^2，x2^2，x1*x2
3rd・・・x1^3，x2^3，x1＊x2^2，x1^2*x2
を考えXTXを取りますと、
x1にはx1^3が交絡、x2にはx2^3が交絡しており、
実は、そのようなモデルは解けないことが分かります。
もともと、中心複合計画は、単峰仮定の計画なのです。
すなわち、ご質問者が解いた式は、
1次モデルそのものという可能性があります。

x1 x2 x1.2 x2.2 x1x2 x1.3 x2.3 x1x2.2 x1.2x2
x1 8.000 0.000 0.000 0.000 0 12.000 0.000 -1.333 0.000
x2 0.000 8.000 0.000 0.000 0 0.000 12.000 0.000 -1.333
x1.2 0.000 0.000 6.666 -1.333 0 0.000 0.000 0.000 0.000
x2.2 0.000 0.000 -1.333 6.666 0 0.000 0.000 0.000 0.000
x1x2 0.000 0.000 0.000 0.000 4 0.000 0.000 0.000 0.000
x1.3 12.000 0.000 0.000 0.000 0 19.999 0.000 -4.000 0.000
x2.3 0.000 12.000 0.000 0.000 0 0.000 19.999 0.000 -4.000
x1x2.2 -1.333 0.000 0.000 0.000 0 -4.000 0.000 2.222 0.000
x1.2x2 0.000 -1.333 0.000 0.000 0 0.000 -4.000 0.000 2.222

ところが、最近のデータサイエンスは、
スパース構造推定ということが可能です。
3次だけのモデルで解きたいのなら、
lasso（ラスー）（罰則付き回帰、L1正則化回帰）という方法で
解かれてはいかがでしょうか。

（上記の行列を求めたRプログラム）
#
x <- data.frame(
matrix(c(
-1, -1,
-1, 1,
1, -1,
1, 1,
-1.4142, 0,
1.4142, 0,
0, -1.4142,
0, 1.4142,
0, 0,
0, 0,
0, 0,
0, 0)
,ncol=2,nrow=12,byrow=T)
)
names(x)=c("x1","x2")
x$x1.2 <- x$x1^2
x$x2.2 <- x$x2^2
x$x1x2 <- x$x1 * x$x2
x <- sweep(x,2,apply(x,2,mean))
round(t(as.matrix(x)) %*% as.matrix(x),digit=3)
x$x1.3 <- x$x1^3
x$x2.3 <- x$x2^3
x$x1x2.2 <- x$x1 * x$x2.2
x$x1.2x2 <- x$x1.2 * x$x2
x <- sweep(x,2,apply(x,2,mean))
round(t(as.matrix(x)) %*% as.matrix(x),digit=3)

yの値を教えて頂ければ、ラスーで行った回帰式を回答します。