統計関連についての覚え書き（※素人なので間違いがあれば教えてください）

DISCLAIMER: No guarantee is given for the correctness of the code herein.

Comparisons among statistical softwares equipped with general mixed model analysis 混合モデルを扱える統計ソフトの比較： SAS(ver. 9.2), SPSS (14.0), R (2.8.1)

split-plotデザインの圃場実験の結果を解析することが最近多い。このデザインでは、block*wholeをrandom effect（変量効果）として解析する必要がある。年次をまたいだ解析の場合は階層がさらに一つ増えsplit-split-plotモデルになり、response = block year whole year*whole sub sub*year sub*whole sub*whole*year がwhole effectに、またrandom effectが、block*year, block*whole+block*year*whole, block*sub + block*sub*year + block*sub*whole + block*sub*year*wholeがとなり、それぞれの階層における残差平方和（error term、誤差項）として用いられる。さらにそれぞれのプロット毎の"個性"を取り除いた解析をしようとすると、そこに共変量（covariate）を入れることになる。また実際には、欠損値（missing value）が生じ、unbalancedなデータセットに対して解析する必要が生じることもある。こういった時に、どの統計ソフトのどういう関数が使えるのか、色々と試行錯誤してみた結果を上の表に示した。

結果は・・・一言で済ますとSASのPROC MIXEDが唯一すべての場合に対応できた。他のソフトだと、欠損値がある場合に重要になってくる、残差平方和の自由度の近似方法が限定されていたり、covariateを主区と副区の誤差項が同じになってしまったり・・・等々、色々と問題が生じた。Rは今後に期待したいが、今のところこの辺りに関する改良の予定は無いようだ。SASは高いだけあってか、SASが！（？）

Covariance analysis in Split-plot design

R

A function "lme" might deal with this problem, but the results differ from those of SAS. I'm still wondering why they differ, but maybe this is because R adopts type I error while SAS uses type III. This conecture is in line with the fact that the results depends on the position of cov (covariate). A single slope is estimated, and tested against sub-plot error term. .

> model<-lme(response ~ cov+whole+sub+whole:sub, random=~1|Block/whole, method="REML", na.action=na.exclude)
> anova (model)
> coefficients(model)

If you use "aov", such like,

> model<-aov(response ~ cov+whole+sub+whole:sub+Error(Block/whole))
> summary (model)
> coefficients(model)

you will get coefficients of the covariate (slope) for every stratum, in this case, Block, whole and sub, hence 3 slopes are estimated. I found that adjusted sums of squares and their DF for whole, error for the whole plot (Block*whole), whole*sub, and error for the subplot (block*whole*sub+block*sub), are identical to the one in a litrature (Federer and King, 2007, Table 11.4 p246), but not for the sub.

SAS

PROC MIXED might deal with this problem, but the results differ from those of lme (R). Change in the position of cov in the model equation does not change the results, probably because this procedure uses type III error by default.

proc mixed data=data method=REML nobound; /*nobound allows variations to go negative, which may be preferred in balanced data*/
class Block whole sub;
model res = Block cov whole sub whole*sub /solution ddfm=CON; /* ddfm defines the way to estimate denomitor degree of freedom. I prefer "containment" (or CON) in balanced case, and "kenwardroger" (or KR) for unbalanced data */
random Block*whole;
LSMEAN whole*sub /CL;
run;

References

Walter T. Federer and Freedom King (2007) Covariance analyses for split plot and split block experiment designs in Variations on split plot and split block experiment designs, pp 239-266.

SASの使い方

Start "cygwin"

startxwin.sh

ssh -X address.of.server

sas_en

Rの使い方

subset: 行列から一部を取り出したいとき

例：subset(data, Week==2) :dataからWeek=2の行のみ抽出する

実践例：

> summary(model<-aov(response ~CO2*Temp*Point+Error(Block/CO2/Temp), subset(data, Week==14)))

> anova(model<-lme(response ~CO2*Temp*Point, random=~1|Block/CO2/Temp,method="ML", na.action=na.exclude, subset(data, Week==14)))

Rで分散分析（ANOVA）

Rでsplit-plotの分析：例えばCO2濃度(CO2)と温度(temperature)が収量(yield)に与える影響を見るため、whole-plot factorがCO2、sub-plot factorがtemperatureで3blockとした場合。

Balanced dataの場合

aov関数を使うと、通常の分散分析を行ってくれる

> data <- read.csv("d:\\....\\....csv", header = T)
> attach(data)
> data
> str(data)
> model <- aov(yields ~ CO2 * Temp + Error(block/CO2))
> summary(model)
> detach(data)

Unbalancedの場合（欠損値missing valueあり）

aovでやってみると、whole-plot, sub-plotそれぞれの誤差に対して主効果を検定してしまうので使えない→他の関数を使う。

• lmeを使った方法（Library: nlme）: Liner Mixed-Effects Model

  > library(nlme)
  > anova(model <- lme(yield ~ CO2*temperature, random= ~1|block/CO2, method = “ML”, na.action=na.exclude))
  > summary(model)

  最尤法を使った解析なので、aovを使った場合と異なる結果になることもある。

• lmerを使った方法（Library: lme4）:Generalized Mixed-Effects Model

> library(lme4)
> lmer(yield ~ CO2 * Temp + (1|Block:CO2), REML = FALSE, family=gaussian, na.action=na.exclude)

lme同様、最尤法を使った解析である。誤差の分布を正規分布以外にしたり、計算プロトコルを変えたりできるらしい。ただしP値を計算してくれないので、検定はできない（おそらく手動で計算すればよいのだが・・・やり方がまだわからず）。

SPSS(Ver. 14.0)でsplit-plotの分散分析

SPSSを使って上記Rと同様のデータを用いて実験結果の分散分析を行ったのでその方法について述べる。また、欠損値が含まれる場合に、うまく誤差の自由度が計算されない、という問題（バグか？）が生じたいのでそれについても報告する。

General Linear Model (GLM)　（一般線型モデル）を用いた解析

変量効果にBlockを、固定効果にCO2とTempを入れて、以下のようにデザインすればよい。通常の分散分析を行ってくれる。以下にSyntax commandを載せる。

UNIANOVA
Response BY CO2 Temp Block
/RANDOM = Block
/METHOD = SSTYPE(3)
/INTERCEPT = INCLUDE
/CRITERIA = ALPHA(.05)
/DESIGN = Block CO2 Block*CO2 Temp CO2*Temp .（Block*CO2はBlock(CO2)と書いても良い）

Mixed model（線型混合モデル）を用いた解析

固定効果にBlock, CO2, Temp, Temp*CO2を、変量効果にBlock*CO2を指定すればよい。ただし、最尤法を使っておりGLMとは異なる結果になることもある。また場合によっては（モデル構造が全く同じでも、"値"によって勝手に判断されてしまう・・・）主区と副区が同じ誤差を持つ構造になってしまう場合があった（この辺は色々議論が分かれるところらしい）。

Syntax:

MIXED
Response BY Block CO2 Temp
/CRITERIA = CIN(95) MXITER(100) MXSTEP(5) SCORING(1)
SINGULAR(0.000000000001) HCONVERGE(0, ABSOLUTE) LCONVERGE(0, ABSOLUTE)
PCONVERGE(0.000001, ABSOLUTE)
/FIXED = Block CO2 Temp CO2*Temp | SSTYPE(3)
/METHOD = REML
/RANDOM Block*CO2 | COVTYPE(VC) .（Block*CO2はBlock(CO2)と書いても良い）

自分で誤差を指定して解析するには？

SPSSは勝手に「最適な」誤差項を探してきて指定してくれるそうだが、自分がやっている限り、階層をもつsplit-plotのような例では、間違った誤差項を指定されてしまう場合もある。そんな時には、ユーザー側で明示的に誤差項を指定する方法としてGLMの中では、TEST Subcommandを使える。

Example:

GLM Response BY CO2 Temp
/TEST = CO2 VS Block*CO2
/TEST = TEMP VS Block*Temp + Block*CO2*Temp
/TEST = CO2*TEMP VS Block*Temp + Block*CO2*Temp

観測値の独立性

分散分析を行う際に最も大事な前提条件に観測値の独立性がある。

サンプル数が増えると誤差の自由度が増え平均平方が小さくなるので有意差が出やすくなる。したがって、独立でないサンプルを独立であると仮定して計算すると、第一種の過誤を起こしやすくなる。分散の均一性、誤差の正規性という他の前提条件に関しては分散分析はそこそこ頑健らしいが、観測値の独立性に関してはシビアになる必要がある。

では観測値が独立であるかどうかは、どう判断すればよいだろうか？

独立である状態とは、「観測値の集合の中の任意の一つの観測値に対して、他の観測値がなんら情報を与えない」状態のことをいうそうだ。といってもわかりにくいので、例を使って説明する。

代表的な「非独立」のサンプルとして、同じ個体から繰り返し測定したデータがある。なぜなら同じ個体からとられたデータは、同じような値をとることが期待されるからである。定義に即した形でこれを記述すると、この場合「ある個体Aでとられた"任意の一つの観測値"に対して、同様に"個体Aでとられた他の観測値"は非常に有力な情報を与えてしまう」ことがわかる。すなわち"ある一つの観測値"は"個体Aでとられた他の観測値"と同じような値をとる、という情報である。したがってこれらの観測値は互いに非独立であると言える。

これを言い換えると、「非独立」の状態とは、「ある観測値の集合があったとき、その中に部分集合ができている状態」と定義できる。「各個体」はそれぞれ一つの集合であるため、一つの個体から何点データをとっても、（おそらくは平均値を用い）一つの値しか独立とはみなせないのである。

こういってもなんだか難しいかもしれないが、考えてみれば高校の時習った、確率の独立性と全く同じことだ。「1」が出ることを期待してサイコロを振ったのに、「5」、「2」、「3」・・・と全然期待した「1」がでない。これだけ振ったのだから、次こそは「1」のハズだ！と思っても、何度サイコロを振ったところで次回「1」が出る確率は「1/6」のままで変わらない。すなわち、これまで振って出た値は、次の試行の結果に対して何ら影響を与えない（別の言い方をすれば、次の試行の結果の予測精度を上げるような情報を与えない）のである。観測値の独立性もこれと同じように考えれば判断がつきやすいのではなかろうか？

例：3反復のブロックを作ったとする。各ブロックの中でいくら観測数を増やしたところで、それを平均した一つの値しか独立とみなせない。なぜなら、そのブロックの中の観測値は互いに似通った値をとると期待されるからである。ただし、観測数を増やすことは、そのブロックの値の代表性を上げるには役に立つ。話はややそれるが、ブロック内の観測数を増やすべきか、あるいはブロックの数自体を増やすべきかは、random effect（変量要因） modelで解析すればヒントが得られる。手短に言えば、ブロック内のばらつきとブロック間のばらつきの程度を比較し、例えばブロック内のバラツキがより大きい場合は、各ブロックないでの観測数を増やして代表性を高めるべき、と考えられる。