読者です 読者をやめる 読者になる 読者になる

R aov関数を分解する

Program R

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

 

1

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

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

 

2

{
    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が追加される。

 

3

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

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

 

4

    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つ以上の値が格納されていれば、処理を中断し、エラーメッセージを表示する。

 

5

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

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

 

6

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

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

 

7

  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を実行する環境をしているもの(詳細要確認)。

 

8

#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がもしも無ければ、エラーを表示する。

 

9

#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関数)の実行後にオプションを保存しておく関数。

 

10

#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に代入する。

 

11

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

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

 

【以下更新予定】

12

    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