R aov関数を分解する

Rでコンソールにaovと入力すると、関数の中身が現れる。普段あまり使わない関数も多数登場しているので、意味を確認することにする。

function (formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, ...)

"formula"は分散分析の対象変数を指定するもの、projectionsは推定を返すかどうか（何の推定？）、qrはQR分解を返すかどうか。contrastsも不明。保留。

{
    Terms <- if (missing(data))　{
　　terms(formula, "Error")} else {
    terms(formula, "Error", data = data)}

missing(data)は、引数dataが省略されていればTRUE、省略されていなければFALSEを返す。

termsはRの様々なデータオブジェクトから"terms"オブジェクトを取り出すための関数。ココによると、Rにはformulaクラスがあり、その中にtermsオブジェクトが格納されている。引数"Error"をつけることで、termsオブジェクトにattr("specials"),attr("specials")$Errorが追加される。

    indError <- attr(Terms, "specials")$Error

上記で作成したTerms(termsオブジェクト)からattr("specials")$Errorを抽出し、indErrorに代入する。("Error"が何を表すのかは、要確認)

    if (length(indError) > 1L) {
        stop(sprintf(ngettext(length(indError), "there are %d Error terms: only 1 is allowed", 
            "there are %d Error terms: only 1 is allowed"), length(indError)), domain = NA)}

もしindErrorに2つ以上の値が格納されていれば、処理を中断し、エラーメッセージを表示する。

    lmcall <- Call <- match.call()
    lmcall[[1L]] <- as.name("lm")
    lmcall$singular.ok <- TRUE

match.call()関数は、RのヘルプのExampleを実行してみると、引数を記録する関数のようだ。引数のうち、最初の引き数(formulaか？)に"lm"という名前を付け、lmcall$sigular.okにTRUEを格納しておく。

    if (projections) {
        qr <- lmcall$qr <- TRUE}
    lmcall$projections <- NULL

projectionがTRUEであれば、lmcall$qrにTUREを格納し、それをqrとしておく。lmcall$projectionにNULLを格納しておく。

  if (is.null(indError)) {
    fit <- eval(lmcall, parent.frame())
    if (projections) 
      fit$projections <- proj(fit)
    class(fit) <- if (inherits(fit, "mlm")) 
      c("maov", "aov", oldClass(fit))
    else c("aov", oldClass(fit))
    fit$call <- Call
    return(fit)
  }　#if節は次に続く

indErrorに何も値が格納されていない場合は、eval関数で文字列(今回ならlmcall)をRで実行する。5でlmcallはformulaが格納されていたから、それを実行して、fitに代入するという意味。parent.frame()はRのヘルプでevalを調べると、evalを実行する環境をしているもの（詳細要確認）。

#if節の続き
  else { #is.null(indError)がFALSEの場合
    if (pmatch("weights", names(match.call()), 0L)) 
      stop("weights are not supported in a multistratum aov() fit")
#if節はまだ続く

indErrorがNULLでない場合は、pmatch関数で引数match.call()のうち、"weights"に該当する引数が、引数match.call()のうち何番目にあるかを調べる。weightsがもしも無ければ、エラーを表示する。

#if節の続き
    opcons <- options("contrasts")
    options(contrasts = c("contr.helmert", "contr.poly"))
    on.exit(options(opcons))
#if節はまだ続く

optionsはグローバルな(リセットしない限り常に有効な)オプションを設定する関数。Rのヘルプによると、contrastsはmodel fittingにおいて使われるオプションのことで、前者がunordered factorsに対して使われるもの、後者がordered factorsに対して使われるもの。"contr.helmert"は第2水準を第1水準と比較し、第3水準を第1水準と第2水準の平均と比較…するオプション。"contr.poly"は"直交多公式系"で比較を行うオプション。(詳細要確認)。on.exitはこの関数(aov関数)の実行後にオプションを保存しておく関数。

#if節の続き 
   allTerms <- Terms
    errorterm <- attr(Terms, "variables")[[1 + indError]]
    eTerm <- deparse(errorterm[[2L]], width.cutoff = 500L, backtick = TRUE)
    intercept <- attr(Terms, "intercept")
    ecall <- lmcall
#if節はまだ続く

Termsは2で作成した通り、formulaのtermsオブジェクト。attr(Terms,"variables")はTermsの変数オブジェクトを表し、その[[1+indError]]番目をerrortermとする、という意味。deparseはwikiによると式を文字列に分解する関数。

atr(Terms,"intercept")はTermsの変数オブジェクト。lmcallはformulaに当たり、それをecallに代入する。

    ecall$formula <- as.formula(paste(deparse(formula[[2L]], width.cutoff = 500L, backtick = TRUE), "~", eTerm, if (!intercept)"- 1"), env = environment(formula))

ecallのformulaオブジェクトに、formula2～eTermを代入する。

【以下更新予定】

    ecall$method <- "qr"
    ecall$qr <- TRUE
    ecall$contrasts <- NULL
    er.fit <- eval(ecall, parent.frame())
    options(opcons)
    nmstrata <- attr(terms(er.fit), "term.labels")
    nmstrata <- sub("^`(.*)`$", "\\1", nmstrata)
    nmstrata <- c("(Intercept)", nmstrata)
    qr.e <- er.fit$qr
    rank.e <- er.fit$rank