7

正規および二重対数プロットのデータに対して、R で線形回帰を実行したいと考えています。

通常のデータの場合、データセットは次のようになります。

lin <- data.frame(x = c(0:6), y = c(0.3, 0.1, 0.9, 3.1, 5, 4.9, 6.2))
plot (lin$x, lin$y)

そこで、データポイント2、3、および4のみの線形回帰の線を引くことを計算したいと思います。

二重対数データの場合、データセットは次のようになります。

data = data.frame(
    x=c(1:15),
    y=c(
        1.000, 0.742, 0.623, 0.550, 0.500, 0.462, 0.433,
        0.051, 0.043, 0.037, 0.032, 0.028, 0.025, 0.022, 0.020
      )
    )
plot (data$x, data$y, log="xy")

ここでは、データセット 1:7 と 8:15 の回帰直線を描きたいと思います。

勾配とy オフセット、およびフィットのパラメーター ( R^2、p-value )を計算できますか?

通常のデータと対数データに対してどのように行われますか?

助けてくれてありがとう、

スヴェン

4

2 に答える 2

13

R では、lm()関数を介して線形最小二乗モデルが適合されます。数式インターフェイスを使用すると、subset引数を使用して、実際のモデルに適合させるために使用されるデータ ポイントを選択できます。次に例を示します。

lin <- data.frame(x = c(0:6), y = c(0.3, 0.1, 0.9, 3.1, 5, 4.9, 6.2))
linm <- lm(y ~ x, data = lin, subset = 2:4)

与える:

R> linm

Call:
lm(formula = y ~ x, data = lin, subset = 2:4)

Coefficients:
(Intercept)            x  
     -1.633        1.500  

R> fitted(linm)
         2          3          4 
-0.1333333  1.3666667  2.8666667

二重ログについては、2 つの選択肢があると思います。i) 上記のように 2 つの別々のモデルを推定するか、ii) ANCOVA を介して推定します。対数変換は、 を使用した式で行われますlog()。

2 つの別々のモデルを介して:

logm1 <- lm(log(y) ~ log(x), data = dat, subset = 1:7)
logm2 <- lm(log(y) ~ log(x), data = dat, subset = 8:15)

または、指標変数が必要な ANCOVA を介して

dat <- transform(dat, ind = factor(1:15 <= 7))
logm3 <- lm(log(y) ~ log(x) * ind, data = dat)

これら 2 つのアプローチが同等かどうかを尋ねるかもしれません。そうですね、モデル係数を介してこれを示すことができます。

R> coef(logm1)
  (Intercept)        log(x) 
-0.0001487042 -0.4305802355 
R> coef(logm2)
(Intercept)      log(x) 
  0.1428293  -1.4966954

したがって、別々のモデルの 2 つの勾配は -0.4306 と -1.4967 です。ANCOVA モデルの係数は次のとおりです。

R> coef(logm3)
   (Intercept)         log(x)        indTRUE log(x):indTRUE 
     0.1428293     -1.4966954     -0.1429780      1.0661152

2つをどのように調整しますか？さて、私がセットアップした方法はind、logm3から推定される値をより直接的に与えるためにパラメータ化されていlogm2ます。と の切片は同じであり、 の係数も同じlogm2です。の係数に相当する値を取得するには、まず切片に対して操作を行う必要があります。logm3log(x)logm1

R> coefs[1] + coefs[3]
  (Intercept) 
-0.0001487042

ここで、 の係数indTRUEは、グループ 1 の平均とグループ 2 の平均の差です。勾配については、次のようになります。

R> coefs[2] + coefs[4]
    log(x) 
-0.4305802

これは、得られたものと同じでありlogm1、グループ 2 の勾配 ( ) に基づいており、coefs[2]グループ 1 の勾配の差 ( ) によって修正されていcoefs[4]ます。

プロットに関しては、簡単な方法はabline()単純なモデルを経由することです。たとえば、通常のデータの例:

plot(y ~ x, data = lin)
abline(linm)

ログ データについては、もう少し工夫が必要かもしれません。ここでの一般的な解決策は、データの範囲を予測し、予測をプロットすることです。

pdat <- with(dat, data.frame(x = seq(from = head(x, 1), to = tail(x,1), 
                                     by = 0.1))
pdat <- transform(pdat, yhat = c(predict(logm1, pdat[1:70,, drop = FALSE]), 
                                 predict(logm2, pdat[71:141,, drop = FALSE])))

累乗することにより、元のスケールでプロットできるyhat

plot(y ~ x, data = dat)
lines(exp(yhat) ~ x, dat = pdat, subset = 1:70, col = "red")
lines(exp(yhat) ~ x, dat = pdat, subset = 71:141, col = "blue")

または対数スケールで：

plot(log(y) ~ log(x), data = dat)
lines(yhat ~ log(x), dat = pdat, subset = 1:70, col = "red")
lines(yhat ~ log(x), dat = pdat, subset = 71:141, col = "blue")

例えば...

この一般的なソリューションは、より複雑な ANCOVA モデルでもうまく機能します。ここで、前と同じように新しい pdat を作成し、インジケーターを追加します

pdat <- with(dat, data.frame(x = seq(from = head(x, 1), to = tail(x,1), 
                                     by = 0.1)[1:140],
                             ind = factor(rep(c(TRUE, FALSE), each = 70))))
pdat <- transform(pdat, yhat = predict(logm3, pdat))

predict()ANCOVA を fit に使用しているため、への 1 回の呼び出しから必要なすべての予測を取得する方法に注目してくださいlogm3。以前と同じようにプロットできます。

plot(y ~ x, data = dat)
lines(exp(yhat) ~ x, dat = pdat, subset = 1:70, col = "red")
lines(exp(yhat) ~ x, dat = pdat, subset = 71:141, col = "blue")
于 2011-06-08T11:07:22.330 に答える
3 
#Split the data into two groups
data1 <- data[1:7, ]
data2 <- data[8:15, ]

#Perform the regression
model1 <- lm(log(y) ~ log(x), data1)
model2 <- lm(log(y) ~ log(x), data2)
summary(model1)
summary(model2)

#Plot it
with(data, plot(x, y, log="xy"))
lines(1:7, exp(predict(model1, data.frame(x = 1:7))))
lines(8:15, exp(predict(model2, data.frame(x = 8:15))))

一般に、データをさまざまなグループに分割し、さまざまなモデルをさまざまなサブセットで実行することは珍しく、おそらく悪い形です。グループ化変数の追加を検討することをお勧めします

data$group <- factor(rep(letters[1:2], times = 7:8))

データセット全体である種のモデルを実行します。たとえば、

model_all <- lm(log(y) ~ log(x) * group, data)
summary(model_all)
于 2011-06-08T11:06:31.957 に答える