目的

重回帰分析を行う。
変数選択を行うには，「sreg 関数」を使用のこと。
（R では lm 関数を使うのがよい）

使用法

mreg(dat, func.name=c("solve", "ginv"))

引数

dat         データ行列。従属変数は最終列におき，それ以外は独立変数と見なす。
func.name   逆行列を計算する関数名（solve か ginv）。
            正規方程式が特異行列になる場合以外はどちらを使っても結果は同じ
            省略時は solve
            ムーア・ペンローズ型一般化逆行列を使うときは ginv

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/mreg.R", encoding="euc-jp")

# 重回帰分析
mreg <- function(    dat,                                    # データ行列
                        func.name=c("solve", "ginv"))           # 逆行列を計算する関数の選択
{
        dat <- subset(dat, complete.cases(dat))                      # 欠損値を持つケースを除く
        n <- nrow(dat)                                               # ケース数
        nc <- ncol(dat)                                              # 列数
        nv <- nc-1                                           # 独立変数の個数（最後の一つは従属変数）
        if (is.null(colnames(dat))) {                           # 変数名が付いていないときには仮の名前を付ける
                colnames(dat) <- paste("Var", 1:nc, sep="")
        }
        r <- cor(dat)                                                # 相関係数行列
        m <- colMeans(dat)[-nc]                                      # 独立変数の平均値ベクトル
        if (match.arg(func.name) == "solve") {
                inverse <- solve
                betas <- inverse(r[-nc, -nc], r[,nc][-nc])   # 標準化偏回帰係数
        }
        else {
                library(MASS)
                inverse <- ginv
                betas <- inverse(r[-nc, -nc]) %*% r[,nc][-nc]        # 標準化偏回帰係数
        }
        variance <- var(dat)*(n-1)                           # 変動共変動行列
        prop <- diag(variance)                                       # 対角成分
        prop <- (prop/prop[nc])[-nc]                         # 偏回帰係数に変換するための係数（独立変数の変動/従属変数の変動）
        b <- betas/sqrt(prop)                                        # 偏回帰係数
        Sr <- variance[,nc][-nc]%*%b                         # 回帰による変動
        St <- variance[nc, nc]                                       # 全変動
        Se <- St-Sr                                          # 誤差変動
        SS <- c(Sr, Se, St)                                  # 平方和（変動）ベクトル
        dfr <- nv                                            # 回帰による変動の自由度
        dfe <- n-nv-1                                                # 誤差変動の自由度
        dft <- n-1                                           # 全変動の自由度
        df <- c(dfr, dfe, dft)                                       # 自由度ベクトル
        Ms <- SS/df                                          # 平均平方ベクトル
        f <- Ms[1]/Ms[2]                                     # F 値
        fvalue <- c(f, NA, NA)
        p <- c(pf(f, dfr, dfe, lower.tail=FALSE), NA, NA)    # P 値
        b0 <- mean(dat[,nc])-sum(b*m)                                # 定数項
        b <- c(b, b0)                                                # 偏回帰係数ベクトル
        inv <- inverse((n-1)*cov(dat)[-nc, -nc])
        SEb <- c(sapply(1:nv, function(i) sqrt(inv[i, i]*Ms[2])), sqrt((1/n+m%*%inv%*%m)*Ms[2]))
        tval <- b/SEb                                                # 偏回帰係数の有意性検定
        pval <- pt(abs(tval), n-nv-1, lower.tail=FALSE)*2    # P 値
        tolerance <- 1/diag(inverse(cor(dat)[-nc, -nc]))     # トレランス
        result <- cbind(b, SEb, tval, pval, c(betas, NA), c(tolerance, NA))
        rownames(result) <- c(colnames(dat)[1:nv], "定数項")
        colnames(result) <- c("偏回帰係数", "標準誤差", "t 値", "P 値", "標準化偏回帰係数", "トレランス")
        R2 <- 1-Se/St                                                # 重相関係数の二乗
        R <- sqrt(R2)                                                # 重相関係数
        R2s <- 1-Ms[2]/Ms[3]                                 # 自由度調整済み重相関係数の二乗
        loglik <- -0.5*n*(log(2*pi)+1-log(n)+log(Se))        # 対数尤度
        AIC <- 2*nc+2-2*loglik                                       # AIC
        Rs <- c("重相関係数"=R, "重相関係数の二乗"=R2, "自由度調整済重相関係数の二乗"=R2s,
                "対数尤度"=loglik, "AIC"=AIC)
        anova <- cbind(SS, df, Ms, fvalue, p)                        # 分散分析表
        rownames(anova) <- c("回帰", "残差", "全体")
        colnames(anova) <- c("平方和", "自由度", "平均平方", "F 値", "P 値")
        return(structure(list(result=result, anova=anova, Rs=Rs), class="mreg"))
}
# print メソッド
print.mreg <- function(      obj,                                    # mreg が返すオブジェクト
                        digits=5)                               # 結果の表示桁数
{
        print(obj$result, digits=digits, na.print="")
        cat("\n回帰の分散分析表\n\n")
        print(obj$anova, digits=digits, na.print="")
        cat("\n")
        sapply(1:length(obj$Rs), function(i) cat(names(obj$Rs)[i], "=", round(obj$Rs[i], digits), "\n"))
}


使用例

> dat <- matrix(c(	# 3 列目が従属変数
+ 	1.2, 1.9, 0.9,
+ 	1.6, 2.7, 1.3,
+ 	3.5, 3.7, 2.0,
+ 	4.0, 3.1, 1.8,
+ 	5.6, 3.5, 2.2,
+ 	5.7, 7.5, 3.5,
+ 	6.7, 1.2, 1.9,
+ 	7.5, 3.7, 2.7,
+ 	8.5, 0.6, 2.1,
+ 	9.7, 5.1, 3.6
+ 	), byrow=TRUE, nc=3)

> (a <- mreg(dat))

       偏回帰係数  標準誤差    t 値       P 値 標準化偏回帰係数 トレランス
Var1      0.20462 0.0075643 27.0506 2.4192e-08          0.67067    0.98358
Var2      0.28663 0.0108015 26.5365 2.7638e-08          0.65793    0.98358
定数項    0.14918 0.0545063  2.7368 2.9050e-02                            

回帰の分散分析表

       平方和 自由度  平均平方   F 値       P 値
回帰 6.671644      2 3.3358218 823.48 4.9319e-09
残差 0.028356      7 0.0040509                  
全体 6.700000      9 0.7444444                  

重相関係数 = 0.99788 
重相関係数の二乗 = 0.99577 
自由度調整済重相関係数の二乗 = 0.99456 
対数尤度 = 15.13807 
AIC = -22.27614 

 解説ページ