Ручное взаимодействие строит линейную регрессию в R

ProgramBox

Ручное взаимодействие строит линейную регрессию в R

Post author:admin
Запись опубликована:4 марта, 2023
Post category:Вопросы по программированию

#r #plot #linear-regression #interaction

#r #график #линейная регрессия #взаимодействие

Вопрос:

Я пытаюсь предсказать среднюю численность животных, замеченных во время разных фаз Луны (фактор), используя логарифмически преобразованные данные о численности (лучше подходят) и некоторые другие переменные. Оказалось, что лучшая модель (наименьший AIC) включает взаимодействие фазы и продолжительности съемки и облачного покрова (оба непрерывные):

 LMoon<-lm(ln~Phase*Duration Clouds, data=abund)

summary(LMoon)

Call:
lm(formula = ln ~ Phase * Duration   Clouds, data = abund)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.75416 -0.46311  0.09522  0.46591  1.85978 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       0.382031   0.876865   0.436 0.664125    
Phase2            2.130065   1.226305   1.737 0.085851 .  
Phase3            1.971060   1.818542   1.084 0.281351    
Phase4            0.608043   1.140122   0.533 0.595146    
Phase5            4.786674   1.151850   4.156 7.44e-05 ***
Phase6            0.958706   1.046831   0.916 0.362238    
Phase7            0.254711   3.425214   0.074 0.940888    
Phase8            0.865995   1.043916   0.830 0.409005    
Duration          0.069153   0.035407   1.953 0.053952 .  
Clouds           -0.004259   0.002401  -1.774 0.079494 .  
Phase2:Duration  -0.087843   0.047818  -1.837 0.069545 .  
Phase3:Duration  -0.089908   0.069652  -1.291 0.200109    
Phase4:Duration  -0.005424   0.046675  -0.116 0.907749    
Phase5:Duration  -0.172016   0.049369  -3.484 0.000768 ***
Phase6:Duration  -0.035597   0.041435  -0.859 0.392583    
Phase7:Duration   0.024084   0.176773   0.136 0.891939    
Phase8:Duration  -0.033424   0.042064  -0.795 0.428963    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7721 on 89 degrees of freedom
Multiple R-squared:  0.3368,    Adjusted R-squared:  0.2176 
F-statistic: 2.825 on 16 and 89 DF,  p-value: 0.0009894

Теперь, из-за этого взаимодействия, мне нужно составить график взаимодействия (CI слишком широк при построении lsmeans).
Я пытался использовать различные функции, упомянутые там, но ни одна из них не работала.
По-видимому, мне нужно рассчитать и построить график вручную, что я и сделал, как это:

 intercepts <- c(coef(LMoon)["(Intercept)"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase2"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase3"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase4"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase5"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase6"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase7"],
                coef(LMoon)["(Intercept)"]   coef(LMoon)["Phase8"])

lines.df <- data.frame(intercepts = intercepts,
                       slopes = c(coef(LMoon)["Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase2:Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase3:Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase4:Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase5:Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase6:Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase7:Duration"],
                                  coef(LMoon)["Duration"] coef(LMoon)["Phase8:Duration"]),
                       Phase2 = levels(abund$Phase))

qplot(x = Duration, y = Sp2, color = Phase, data = abund)   
  geom_abline(aes(intercept = intercepts, 
                  slope = slopes, 
                  color = Phase), data = lines.df)

График, который я получаю, неверен, поскольку значения y находятся в исходном истинном масштабе, но строки основаны на lm, который использует данные, преобразованные в журнал.

количество графиков взаимодействия, продолжительность, фазы Луны

Для обратного преобразования кто-то сказал мне, что в итоге я не получу прямых линий. Вместо того, чтобы использовать abline (), я должен создать набор, например, из 100 новых значений x, которые охватывают диапазон данных о продолжительности, и использовать коэффициенты для вычисления ваших прогнозируемых значений y. Затем постройте их с помощью lines(), и это должно выглядеть как плавная кривая.

И вот тут я заблудился.

Итак, я создал набор новых значений x для диапазона продолжительности опроса (мин. 15 максимум 45 мин): dur2 <- seq(from = 15, to = 45, length.out=100)

Затем, как только я получу эти значения, я должен получить прогнозируемое значение y для каждого значения x, используя коэффициенты моего LM. После этого выполните обратное преобразование значений y в исходный масштаб. А затем, используя значения x и преобразованные обратно значения y, добавьте линии к графику.

Как мне теперь точно получить прогнозируемые значения? Я не могу использовать какой-либо тип / функцию pred, я все это перепробовал. Это просто не работает с моей моделью, поэтому единственный способ — вручную, но я понятия не имею, как это сделать…

Надеюсь, что кто-нибудь сможет мне в этом помочь, я пытался неделями и к настоящему времени в отчаянии, близок к тому, чтобы сдаться.

Приветствия!

PS Здесь данные:

 > dput(subset(abund, Phase %in% c("Phase1", "Phase2")))

structure(list(Year = integer(0), Date = structure(integer(0), .Label = c("01/08/2009", 
"01/08/2016", "02/07/2019", "02/08/2009", "02/08/2012", "02/08/2016", 
"02/09/2007", "03/08/2007", "03/08/2009", "03/08/2014", "03/08/2015", 
"04/07/2019", "04/08/2009", "04/08/2013", "05/08/2009", "05/08/2014", 
"05/08/2015", "06/07/2008", "06/07/2019", "07/08/2009", "08/07/2010", 
"09/07/2010", "09/08/2015", "10/08/2009", "11/08/2009", "12/08/2009", 
"13/08/2009", "13/08/2014", "14/08/2009", "14/08/2012", "16/07/2006", 
"18/07/2009", "18/08/2015", "19/07/2011", "20/08/2009", "21/07/2011", 
"21/09/2009", "22/07/2011", "22/07/2016", "22/07/2017", "23/07/2007", 
"23/07/2016", "23/07/2017", "24/07/2017", "25/07/2007", "25/07/2010", 
"25/07/2017", "25/08/2016", "26/07/2010", "26/07/2011", "27/07/2006", 
"27/07/2011", "27/07/2012", "28/07/2016", "29/06/2019", "29/07/2005", 
"29/07/2009", "29/07/2010", "29/07/2016", "29/07/2019", "30/07/2005", 
"30/07/2007", "30/07/2016", "30/08/2005", "31/07/2005", "31/07/2009", 
"31/07/2014", "31/07/2016"), class = "factor"), NrSurvey = integer(0), 
    Duration = integer(0), Sp2 = integer(0), Phase = structure(integer(0), .Label = c("1", 
    "2", "3", "4", "5", "6", "7", "8"), class = "factor"), Clouds = integer(0), 
    Visibility = integer(0), ln = numeric(0)), row.names = integer(0), class = "data.frame")

1. tl; dr но смотрите: nature.com/articles/npre.2010.4136.1 почему вы не должны преобразовывать данные подсчета в журнал

2. Это действительно легко сделать с помощью пакета emmeans.

3. @Dylan_Gomes это правда, но у меня нет нулей в моем наборе данных, поскольку я фокусируюсь только на положительных опросах, т.Е. Там, Где были замечены животные. Я просто хочу посмотреть, отличается ли количество животных в зависимости от фаз Луны. Для отсутствия присутствия я использую разные модели

4. @RussLenth не могли бы вы уточнить? Я уже пробовал emmeans раньше, однако CI слишком широк из-за взаимодействия

5. @Cathrin, избавление от нулей также приводит к потере информации, поскольку отсутствие — это информация. Возможно, стоит взглянуть на пуассоновскую или отрицательную биномиальную модель, а не на то, что вы делаете. Но это только мое мнение.

Ответ №1:

Используется predict для получения прогнозируемых значений. Не вычисляйте вручную.

Используется expand.grid() для создания фрейма данных всех комбинаций вашей dur2 последовательности и других предикторов со значением (значениями), которые вы хотите отобразить. Что-то вроде этого:

 prediction_data = expand.grid(
  Duration = dur2,
  Phase= unique(abund$Phase),
  Clouds = mean(abund$Clouds) # Hold Clouds constant at some value
)

# column names in prediction_data need to match those in the model formula
prediction_data$pred = predict(LMoon, newdata = prediction_data)
prediction_data$pred_orig = exp(prediction_data$pred)

# plot
ggplot(prediction_data, aes(x = Duration, y = pred_orig, color = Phase))  
  geom_line()  
  geom_point(data = abund)

Что-то подобное должно сработать.

Еще один хороший вариант — использовать broom::augment для генерации прогнозов. Это может легко дать стандартные ошибки и остатки для каждой точки прогнозирования.

 library(broom)
prediction_data = augment(LMoon, newdata = prediction_data)

1. Убедитесь prediction_data , что у столбца есть имя Phase , и все должно быть в порядке. Если вы опубликуете воспроизводимую выборку данных (возможно dput(subset(abund, Phase %in% c("Phase1", "Phase2"))) , или, возможно, меньшее подмножество, тогда мы действительно сможем запустить код на чем-то, похожем на ваши данные, что очень полезно для отладки.

2. Большое спасибо, это сработало!! — Кроме построения графика, он продолжает выдавать мне ошибку FUN(X[[i]], ...) : object 'pred_orig' not found . Однако при запуске только строки ggplot(prediction_data, aes(x = Duration, y = pred_orig)) без geom_line / point() ошибка не возникает. Есть идеи?

3. Также заметил, что график пуст. Это не при использовании base = 10 в log() Перед запуском модели и т. Д., Тогда график дает мне, по крайней мере, несколько точек и линий, но все равно не с geom_line / point()

4. Вам нужно будет создать pred_orig столбец, как в моем ответе. Понятия не имею о пустом графике, но если вы поделитесь некоторыми примерами данных, я мог бы попытаться взглянуть. Смотрите мой предыдущий комментарий о том, как обмениваться образцами данных.

5. Привет, у меня возникли другие проблемы с графиком взаимодействия — с двоичными данными присутствия / отсутствия. График соединяет каждую точку между верхним и нижним пределом, есть идеи, почему это так? Я использовал ваше предложение выше, но это было для подсчета данных. Работа с glm()