Как создать диаграмму рассеяния, отображающую несколько предсказаний в линейной регрессии трехсторонних взаимодействий с использованием базового R - PullRequest
Новая
       47

Как создать диаграмму рассеяния, отображающую несколько предсказаний в линейной регрессии трехсторонних взаимодействий с использованием базового R

0 голосов
/ 09 октября 2019

Получив модель накопления биомассы во времени в трехстороннем взаимодействии, включающем факториальную переменную с 2 уровнями (традиционным и интенсивным) и непрерывную переменную (EC), я попытался представить результат на рисунке, используя диаграмму рассеянияэто будет отображать линии регрессии для четырех разных случаев: традиционное лечение и высокая электропроводность;традиционное лечение и низкая электропроводность;интенсивная обработка и высокая электропроводность;интенсивное лечение и низкая электропроводность.

Я пытался использовать ggplot, но он не обеспечивает достаточной гибкости в том, как я представляю свой результат. Я хочу использовать базовую букву R, чтобы точно представить диапазоны электропроводности, которые я хочу использовать, и более точно персонализировать цветовую палитру и условные обозначения.

Для этого я создал 4 набора данных для построения регрессии. линии из модели, примененной к набору данных 

data
    age  logeage    use     EC         AGB
1   9.5 2.251292 intens 0.0825   0.0000000
2  65.0 4.174387   trad 0.2625 237.6903265
3  11.0 2.397895 intens 0.0850   0.4857754
4  26.0 3.258097   trad 0.2350 263.9812024
5  29.0 3.367296   trad 0.4300 106.6392453
6  13.5 2.602690 intens 0.5100   0.0000000
7  29.0 3.367296   trad 0.4400 241.2897229
8  13.5 2.602690 intens 0.2800   0.0000000
9  12.5 2.525729 intens 0.2000   1.9553807
10 21.0 3.044522 intens 0.1300   5.4259407
11  1.0 0.000000   trad 0.2600   0.0000000
13 25.0 3.218876   trad 0.2900  98.0180916
14 30.0 3.401197   trad 0.7600 197.2887138
15 31.0 3.433987   trad 0.3900 254.8277552
16 66.0 4.189655   trad 0.2400 135.1967346
17 30.5 3.417727   trad 0.2100 217.6831019
18 12.0 2.484907 intens 0.0700   0.0000000
19 48.0 3.871201   trad 0.3400 105.1035670
20 15.0 2.708050 intens 0.2700  12.8276905
21 30.5 3.417727   trad 0.4400 167.0891661
22 32.0 3.465736 intens 0.5800 108.9973034
23 71.0 4.262680   trad 0.3700 135.9036366
24 35.5 3.569533   trad 0.1300  97.1231771
25 30.0 3.401197   trad 0.3100 107.0933743
27 35.0 3.555348   trad 0.1500 146.8864254
28 69.0 4.234107   trad 0.1900 267.8302222
29 12.0 2.484907 intens 0.0600  21.4084421
34 30.5 3.417727   trad 0.1300  72.9122228
35 19.5 2.970414   trad 0.1400  17.3043656
37 11.5 2.442347   trad 0.2800  91.5393234
38 28.5 3.349904 intens 0.2300 198.6790722
39 47.0 3.850148   trad 0.1100 222.4857276
40 46.0 3.828641   trad 0.1900 171.3323051
41 43.0 3.761200   trad 0.1100 184.9173687
42  9.0 2.197225 intens 0.2900  10.7360386
43 22.0 3.091042   trad 0.3700 149.7812194
49 28.5 3.349904   trad 0.2600 112.7302377
50 51.0 3.931826 intens 0.2200  21.9547804
55 19.0 2.944439   trad 0.1200  58.6177275
56 28.0 3.332205 intens 0.0700  31.1173969
57 47.0 3.850148   trad 0.1100 194.4897458
58 45.0 3.806662   trad 0.1100  90.0451401
59  7.5 2.014903 intens 0.0900   2.4621935
60 10.0 2.302585 intens 0.0400  12.2489902
61 36.0 3.583519   trad 0.0500  95.4059899

rbPal <- colorRampPalette(c('red','yellow'))
data$colour<- rbPal(10)[as.numeric(cut(data$EC, breaks = 10))]
data$use<- trim_levels(data$use)

  par(mfrow = c(1,1), mar= c(2.5,2.5,1,0.5), yaxt= "s",xaxt= "s", las = 0, cex = 1)
  plot(AGB~age, data = data, pch= c(1,8)[as.numeric(use)], cex= 1, col= colour,
       xlab="Time", ylab="biomass", ylim = c(0,median(pristine$AGB)+10), mgp=c(1.5,0.4,0))

  # prediction traditional treatment
  mdl.r <- lm(AGB ~ logeage*use*EC, data = data)
  data.low<- data.frame(logeage=seq(from= min(data$logeage[data$EC < mean(
    data$EC)]),to=max(data$logeage[data$EC < mean(data$EC)]),length=50), EC = seq(
      from= min(data$EC), to=mean(data$EC), length=50), use = "trad")
  data.high<- data.frame(logeage=seq(from= min(data$logeage[data$EC > mean(
    data$EC)]), to=max(data$logeage[data$EC > mean(data$EC)]), length=50), EC = seq(
      from= mean(data$EC), to=max(data$EC), length=50), use = "trad")

  pred.H<-  predict(mdl.r, newdata=data.high, type="response", se = T)
  pred.L<-  predict(mdl.r, newdata=data.low, type="response", se = T)
  lines(exp(data.high$logeage), pred.H$fit, lty=3, col="red", lwd = 5)
  lines(exp(data.low$logeage), pred.L$fit, lty=3, col="orange", lwd = 3)

  # prediction intensive treatment
  mdl.r <- lm(AGB ~ logeage*use*EC, data = data)
  data.low<- data.frame(logeage=seq(from= min(data$logeage[data$EC < mean(
    data$EC)]),to=max(data$logeage[data$EC < mean(data$EC)]),length=50), EC = seq(
      from= min(data$EC), to=mean(data$EC), length=50), use = "intens")
  data.high<- data.frame(logeage=seq(from= min(data$logeage[data$EC > mean(
    data$EC)]), to=max(data$logeage[data$EC > mean(data$EC)]), length=50), EC = seq(
      from= mean(data$EC), to=max(data$EC), length=50), use = "intens")
  pred.H<-  predict(mdl.r, newdata=data.high, type="response", se = T)
  pred.L<-  predict(mdl.r, newdata=data.low, type="response", se = T)
  lines(exp(data.high$logeage), pred.H$fit, lty=1, col="red", lwd = 5)
  lines(exp(data.low$logeage), pred.L$fit, lty=1, col="orange", lwd = 3)

Из 4 прогнозов, 3, по-видимому, соответствуют ожидаемому результату: логарифмическое увеличение AGB во времени, однако прогноз для традиционного лечения + высокий ECпоказывает совершенно неожиданную тенденцию, показывая артефакт, который может возникнуть из-за ошибки в наборе данных, созданном для прогноза. Я ожидаю, что все кривые вырастут в геометрической прогрессии, как и три другие.

...