Возникли проблемы с сохранением подобранных и наблюдаемых значений в модели lm?

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

Я добился сохранения некоторых из них, включая коэффициенты R2 и Pvalue. Однако, когда я сохраняю подобранные и наблюдаемые значения в data.frame с именем obsFit, он возвращает подобранные и наблюдаемые значения только одного категориального уровня, поскольку мне нужны уровни 7.

Есть идеи и мысли?

Извините за жесткое кодирование R, я все еще новичок, который находится в процессе обучения.

Вот патрон code, который я использовал.

# Wilayah 1 vs MADA4 model
> Y1 <- yield[ ,2:5]
> mada4 <- stns[["S1"]][["Mada4_sum"]]
> unique(mada4$aggre_per)
[1] "3"     "4"     "5"     "6"     "7"     "mamjj"
[7] "amj"  
> 
> # create an empty matrix to fill results
> df <- matrix(rep(NA, 9))
> 
> for (agg.per in unique(mada4$aggre_per)) {
+   
+   subdata <- subset(mada4, aggre_per == agg.per)
+   subdata <- cbind(subdata, Y1$I_S1)
+   names(subdata)[4] <- "I_S1"
+   
+   #save model summary
+   mod <- lm('I_S1 ~ value', subdata)
+   slope     <- summary(mod)[[4]][2,1]
+   intercept <- summary(mod)[[4]][1,1]
+   r2        <- summary(mod)[[8]]
+   adjr2     <- summary(mod)[[9]]
+   fstat     <- summary(mod)[[10]][[1]]
+   pval      <- summary(mod)[[4]][2,4]
+   pearson   <- cor(subdata$I_S1,subdata$value )
+   res <- c("Mada4_sum", agg.per, slope, intercept, r2, adjr2, fstat, pval, pearson)
+   df <- cbind(df, res)
+   
+   #save the fitted and observed yields as data.frame
+   obsFit <- data.frame(cbind("Mada4_sum", agg.per,subdata$I_S1, fitted(mod)))
+   
+ }
> 
> df <- data.frame(t(df[,-1]))
> df[,-2:-1] <- apply (df[,-2:-1], 2, as.numeric)
> df[,-2:-1] <- apply (df[,-2:-1], 2, round, digit =2)
> rownames(df) <- 1:nrow(df)
> names(df) <-c( "WI",
+                "aggre_per",
+                "slope",
+                "intercept",
+                "R2",
+                "adjustedR2",
+                "F-stat",
+                "p-value",
+                "correlation")
> head(df)
         WI aggre_per slope intercept   R2 adjustedR2
1 Mada4_sum         3 -1.17   6230.91 0.09       0.03
2 Mada4_sum         4 -0.27   6128.81 0.00      -0.07
3 Mada4_sum         5 -1.04   6240.46 0.02      -0.05
4 Mada4_sum         6 -2.88   6507.82 0.14       0.08
5 Mada4_sum         7 -1.12   6296.72 0.05      -0.02
6 Mada4_sum     mamjj -0.89   6760.92 0.18       0.12
  F-stat p-value correlation
1   1.40    0.26       -0.30
2   0.07    0.80       -0.07
3   0.23    0.64       -0.13
4   2.26    0.15       -0.37
5   0.76    0.40       -0.23
6   3.12    0.10       -0.43
> #check levels of `aggre_per`
> unique(df$aggre_per)
[1] "3"     "4"     "5"     "6"     "7"     "mamjj"
[7] "amj"  
> 
> names(obsFit) <- c("WI", "Aggr.per", "Observed", "Fitted")
> print(obsFit)
           WI Aggr.per Observed           Fitted
103 Mada4_sum      amj     5204 6019.05130014577
104 Mada4_sum      amj     5824 6061.66361454587
105 Mada4_sum      amj     5481 6129.49546195825
106 Mada4_sum      amj     5587 5962.52476063545
107 Mada4_sum      amj     6260 6124.27762754192
108 Mada4_sum      amj     5771   6017.312022007
109 Mada4_sum      amj     5899 5888.60543973734
110 Mada4_sum      amj     6338 6189.50055774613
111 Mada4_sum      amj     6340 6023.39949549272
112 Mada4_sum      amj     6013 6146.01860427665
113 Mada4_sum      amj     6075 6221.67720331355
114 Mada4_sum      amj     6484 5959.91584342728
115 Mada4_sum      amj     6771 5998.17996248043
116 Mada4_sum      amj     6493 6249.50565353401
117 Mada4_sum      amj     6386 6290.37868979533
118 Mada4_sum      amj     6465 6109.49376336229
> 

Вот dput моих данных?

> dput(mada4)
structure(list(year = c(2001, 2002, 2003, 2004, 2005, 2006, 2007, 
2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2001, 2002, 
2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 
2014, 2015, 2016, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 
2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2001, 2002, 2003, 
2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 
2015, 2016, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 
2010, 2011, 2012, 2013, 2014, 2015, 2016, 2001, 2002, 2003, 2004, 
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 
2016, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 
2011, 2012, 2013, 2014, 2015, 2016), aggre_per = c("3", "3", 
"3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", "3", 
"3", "4", "4", "4", "4", "4", "4", "4", "4", "4", "4", "4", "4", 
"4", "4", "4", "4", "5", "5", "5", "5", "5", "5", "5", "5", "5", 
"5", "5", "5", "5", "5", "5", "5", "6", "6", "6", "6", "6", "6", 
"6", "6", "6", "6", "6", "6", "6", "6", "6", "6", "7", "7", "7", 
"7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", "7", 
"mamjj", "mamjj", "mamjj", "mamjj", "mamjj", "mamjj", "mamjj", 
"mamjj", "mamjj", "mamjj", "mamjj", "mamjj", "mamjj", "mamjj", 
"mamjj", "mamjj", "amj", "amj", "amj", "amj", "amj", "amj", "amj", 
"amj", "amj", "amj", "amj", "amj", "amj", "amj", "amj", "amj"
), value = c(189, 6, 212, 81, 101, 84, 163, 198, 143, 128, 456, 
93, 37, 23, 58, 0, 210, 196, 70, 122, 132, 184, 387, 57, 232, 
153, 11, 325, 228, 75, 41, 31, 141, 95, 219, 217, 178, 154, 139, 
112, 164, 62, 106, 202, 101, 123, 108, 252, 175, 186, 110, 252, 
95, 190, 150, 161, 125, 165, 176, 67, 221, 63, 65, 139, 100, 
267, 207, 202, 48, 322, 240, 250, 270, 298, 159, 164, 178, 29, 
94, 164, 815, 750, 818, 874, 554, 934, 1079, 778, 934, 806, 908, 
851, 765, 313, 366, 586, 526, 477, 399, 591, 405, 528, 676, 330, 
521, 380, 293, 594, 550, 261, 214, 422)), row.names = c(1L, 2L, 
3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 
18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 
31L, 32L, 33L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 
45L, 46L, 47L, 48L, 49L, 50L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 
59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 69L, 70L, 71L, 72L, 
73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 86L, 
87L, 88L, 89L, 90L, 91L, 92L, 93L, 94L, 95L, 96L, 97L, 98L, 99L, 
100L, 101L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L, 111L, 
112L, 113L, 114L, 115L, 116L, 117L, 118L), class = "data.frame")
> 
> dput(yield)
structure(list(date = structure(c(11323, 11688, 12053, 12418, 
12784, 13149, 13514, 13879, 14245, 14610, 14975, 15340, 15706, 
16071, 16436, 16801), class = "Date"), I_S1 = c(5204, 5824, 5481, 
5587, 6260, 5771, 5899, 6338, 6340, 6013, 6075, 6484, 6771, 6493, 
6386, 6465), II_S1 = c(5081, 5674, 5589, 5809, 6055, 6231, 5728, 
5971, 6268, 5999, 6015, 6624, 6827, 6453, 6382, 6340), III_S1 = c(5071, 
5349, 5703, 5451, 5764, 5528, 5629, 5395, 5797, 5785, 5887, 5581, 
5899, 6179, 6086, 5734), IV_S1 = c(5515, 5491, 5969, 5971, 6265, 
6242, 6153, 6188, 6311, 6447, 6213, 6410, 6744, 6754, 6247, 6415
), I_S2 = c(4735, 5227, 5415, 5894, 3294, 5113, 6056, 6216, 5752, 
5758, 5540, 5418, 6461, 5346, 6437, 5633), II_S2 = c(4837, 5209, 
5111, 5482, 3203, 5389, 5209, 6041, 5617, 5927, 4993, 5342, 6380, 
5122, 6320, 5470), III_S2 = c(4720, 4986, 4846, 4720, 3968, 4912, 
4437, 5675, 5368, 5638, 4892, 4849, 6117, 4996, 5586, 4666), 
    IV_S2 = c(5054, 5584, 5546, 5666, 5116, 6288, 6089, 6529, 
    5880, 6156, 5263, 5358, 6804, 5212, 6174, 5560)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -16L))

Вы перезаписываете obsFit во время каждой итерации цикла. Если вы хотите сохранить все выходные данные, вам нужно будет поместить их в список фреймов данных.

Allan Cameron 09.04.2022 12:48

Пожалуйста, также помогите, если вы можете написать мой код короче, так как я чувствую, что я жестко запрограммирован. Спасибо

Mukhtar 09.04.2022 13:02
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
1
2
33
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Вы можете сделать это в by, который разбивает данные по условию и применяет функцию, а затем rbind результирующий список. obsFit, через которое вы могли бы пройти как attribute.

res <- by(mada4, mada4$aggre_per, \(x) {
  x <- cbind(x, yield[ ,2:5])
  sm <- summary(mod <- lm(I_S1 ~ value, x))
  o <- setNames(as.vector(sm$coefficients[, c(1, 4)])[c(1, 2, 4)],
                c('Intercept', 'slope', 'pval'))
  o <- c(o, sapply(sm[c('r.squared', 'adj.r.squared', 'fstatistic')], `[`, 1))
  o <- c(o, pearson=with(x, cor(I_S1, value)))[c(2, 1, 4, 5, 6, 3, 7)]
  ft <- cbind(x, x$I_S1, fitted(mod))
  return(`attr<-`(o, 'obsFit', ft))
})

Полученные результаты

do.call(rbind, res) |> signif (4)
#         slope Intercept r.squared adj.r.squared fstatistic.value    pval  pearson
# 3     -1.1680      6231   0.09105      0.026130          1.40200 0.25600 -0.30170
# 4     -0.2730      6129   0.00464     -0.066460          0.06526 0.80210 -0.06812
# 5     -1.0350      6240   0.01636     -0.053900          0.23290 0.63680 -0.12790
# 6     -2.8780      6508   0.13910      0.077660          2.26300 0.15470 -0.37300
# 7     -1.1220      6297   0.05134     -0.016420          0.75760 0.39880 -0.22660
# amj   -0.8696      6476   0.07058      0.004188          1.06300 0.32000 -0.26570
# mamjj -0.8889      6761   0.18210      0.123700          3.11700 0.09926 -0.42670

obsFits

obsFits <- lapply(res, attr, 'obsFit')
obsFits$`3`  ## first aggre_per `3`
#    year aggre_per value I_S1 II_S1 III_S1 IV_S1 x$I_S1 fitted(mod)
# 1  2001         3   189 5204  5081   5071  5515   5204    6010.131
# 2  2002         3     6 5824  5674   5349  5491   5824    6223.904
# 3  2003         3   212 5481  5589   5703  5969   5481    5983.264
# 4  2004         3    81 5587  5809   5451  5971   5587    6136.292
# 5  2005         3   101 6260  6055   5764  6265   6260    6112.929
# 6  2006         3    84 5771  6231   5528  6242   5771    6132.788
# 7  2007         3   163 5899  5728   5629  6153   5899    6040.503
# 8  2008         3   198 6338  5971   5395  6188   6338    5999.618
# 9  2009         3   143 6340  6268   5797  6311   6340    6063.866
# 10 2010         3   128 6013  5999   5785  6447   6013    6081.389
# 11 2011         3   456 6075  6015   5887  6213   6075    5698.234
# 12 2012         3    93 6484  6624   5581  6410   6484    6122.274
# 13 2013         3    37 6771  6827   5899  6744   6771    6187.691
# 14 2014         3    23 6493  6453   6179  6754   6493    6204.045
# 15 2015         3    58 6386  6382   6086  6247   6386    6163.160
# 16 2016         3     0 6465  6340   5734  6415   6465    6230.913

Примечание: R >= 4.1 используется.

Другие вопросы по теме