目的

Excel にある二変量統計関数を R で定義する

使用法

correl(x, y)
covar(x, y)
forecast(data, y, x)
growth(y, x, data, one=FALSE)
intercept(y, x)
logest(y, x, one=FALSE)
pearson(x, y)
rsq(y, x)
slope(y, x, zero=FALSE)
steyx(y, x)
trend(y, x, data, zero=FALSE)

引数

x は独立変数
y は従属変数
data は予測値を求める独立変数の値
growth，logest 関数の引数の one は FALSE または TRUE の値を取る。y = a bx において，TRUE の場合には a = 1 として b のみを求める。FALSE の場合は a，b ともに求める。FALSE の場合には one=FALSE は省略できる。
slope，trend 関数の引数の zero は FALSE または TRUE の値を取る。y = a+b x において，TRUE の場合には a = 0 として b のみを求める（原点を通る回帰直線）。FALSE の場合は a，b ともに求める。FALSE の場合には zero=FALSE は省略できる。

ソース

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

# Excel の二変量統計関数
correl <- cor                                                # 単に名前の違い。ただし，cor(x, y) の場合のみ

covar <- function(x, y)                                      # 共分散（不偏共分散ではない）
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        n <- length(x)                                       # ケース数
        var(x, y)*(n-1)/n                               # Excel は共変動をデータ数（n）で割ったものとして定義している
}

forecast <- function(data, y, x)                     # R の特徴で，data はスカラーでもベクトルでもかまわない
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        mean(y)+var(x, y)/var(x)*(data-mean(x))         # y，x について回帰直線を求め，独立変数の値が data のときの予測値を求める
}

growth <- function(y, x, data, one = FALSE)          # R の特徴で，data はスカラーでもベクトルでもかまわない
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        stopifnot(all(y > 0))                           # 従属変数はすべて正の値を取らなければならない
        y <- log(y)                                  # y = a*b^x の，両辺の対数をとり，直線回帰に持ち込む
        if (one) {                                      # y，x について y = b^x にあてはめ，独立変数の値が data のときの予測値を求める
                b <- sum(x*y)/sum(x^2)                       # 原点を通る直線
                const <- 0
        }
        else {                                          # y，x について y = a*b^x にあてはめ，独立変数の値が data のときの予測値を求める
                b <- var(x, y)/var(x)                        # 切片を持つ直線
                const <- mean(y)-b*mean(x)
        }
        exp(const+b*data)                               # 予測値を計算し，その指数をとって元の尺度に戻す
}

intercept <- function(y, x)
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        mean(y)-var(x, y)/var(x)*mean(x)                # 回帰直線の切片
}

logest <- function(y, x, one = FALSE)
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        stopifnot(all(y > 0))                           # 従属変数はすべて正の値を取らなければならない
        y <- log(y)                                  # y = a*b^x の，両辺の対数をとり，直線回帰に持ち込む
        if (one) {                                      # y，x について y = b^x にあてはめ，b を求める
                b <- sum(x*y)/sum(x^2)                       # 原点を通る直線
                const <- 0
        }
        else {                                          # y，x について y = a*b^x にあてはめ，a と bを求める
                b <- var(x, y)/var(x)                        # 切片を持つ直線
                const <- mean(y)-b*mean(x)
        }
        list(model="a*b^x", result=c(a=exp(const), b=exp(b)))
}

pearson <- cor                                               # 単に名前の違い

rsq <- function(y, x)
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        cor(y, x)^2                                     # 単相関係数の二乗
}

slope <- function(y, x, zero = FALSE)
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        ifelse(zero, sum(x*y)/sum(x^2),                 # zero=TRUE  のときには，原点を通る回帰直線の傾きを計算する
                var(x, y)/var(x))                       # zero=FALSE のときには，切片を持つ回帰直線の傾きを計算する
}

steyx <- function(y, x)
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        n <- length(x)
        sqrt((var(y)-var(x, y)^2/var(x))*(n-1)/(n-2))   # 回帰直線において，回帰からの標本標準偏差（残差に対する平均平方の平方根）を計算する
}

trend <- function(y, x, data, zero = FALSE)          # R の特徴で，data はスカラーでもベクトルでもかまわない
{
        OK <- complete.cases(x, y)                   # 二変数ともに欠損値を持たないケースを選択する
        x <- x[OK]
        y <- y[OK]
        
        ifelse(zero, sum(x*y)/sum(x^2)*data,            # zero=TRUE  のときには，原点を通る回帰直線による予測値を計算する
                intercept(y, x)+var(x, y)/var(x)*data)  # zero=FALSE のときには，切片を持つ回帰直線による予測値を計算する
}


使用データ例

x <- c(3, 5, 4, 8, 2, 1)	# 独立変数
y <- c(2, 4, 5, 7, 1, 6)	# 従属変数

使用例と出力結果例

> correl(x, y)
[1] 0.4925158

> covar(x, y)
[1] 2.361111

> forecast(6, y, x)
[1] 5.162162

> growth(y, x, 6)
[1] 4.677107

> growth(y, x, 6, one=TRUE)
[1] 5.228738

> intercept(y, x)
[1] 2.405405

> logest(y, x)
$model
[1] "a*b^x"

$result
       a        b 
2.010295 1.151117 

> logest(y, x, one=TRUE)
$model
[1] "a*b^x"

$result
       a        b 
1.000000 1.317446 

> pearson(x, y)
[1] 0.4925158

> rsq(y, x)
[1] 0.2425718

> slope(y, x)
[1] 0.4594595

> slope(y, x, zero=TRUE)
[1] 0.9243697

> steyx(y, x)
[1] 2.254125

> trend(y, x, 6)
[1] 5.162162

> trend(y, x, 6, zero=TRUE)
[1] 5.546218

 解説ページ