Линейная регрессия в R (нормальные и логарифмические данные)

Question

Я хочу провести линейную регрессию в 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-значение )

Как это делается для нормальных и логарифмических данных?

Спасибо за помощь,

Sven

Answer 1

В 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

Что касаетсядвойной журнал, у вас есть два варианта, я думаю;i) оценить две отдельные модели, как мы делали выше, или ii) оценить через ANCOVA.Преобразование журнала выполняется в формуле с использованием log().

Через две отдельные модели:

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)

Вы можете спросить, эквивалентны ли эти два подхода?Ну, они есть, и мы можем показать это через коэффициенты модели.

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

Таким образом, два уклона -0,4306 и -1,4967 для отдельных моделей.Коэффициенты для модели ANCOVA:

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

Как мы согласовываем эти два?Хорошо, способ, которым я настроил ind, logm3, параметризован, чтобы дать более прямые значения, оцененные из logm2;точки пересечения logm2 и logm3 такие же, как и коэффициенты для log(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, чтобы соответствовать 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")

Answer 2

#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)