2
\$\begingroup\$

PROBLEM STATEMENT
I have a code to treat agronomic data. I have a data.frame called yieldwith different information (dry.weight / infected pod / healthy pod). These data are collected in the field every month. I would like to have the total yield per year and per tree (plotnb).

Another important thing I'm seeking is the total number of pods produced (either diseased or healthy) the problem I have is that some pods are preleved to make some disease resistance test and therefore are not counted as either healthy or diseased.

It is a routine code and (with whole data) it takes around 1.5 minutes to calculate and I find it too much for what I'm calculating so I'm sure there is a cleverest way to do it.

I'm seeking a code review here to see what you guys think of my code and if there is a way to optimize it and I will take any advice you could give me to improve it!

REPRODUCIBLE DATA

yield <- structure(list(plotnb = c(49L, 49L, 49L, 49L, 49L, 49L, 49L, 
                      49L, 49L, 49L, 49L, 89L, 158L, 158L, 158L, 158L, 158L, 158L, 
                      158L, 158L, 158L, 159L, 159L, 249L, 249L, 249L, 318L, 318L, 318L, 
                      326L, 326L, 326L, 326L, 326L, 349L, 349L, 408L, 421L, 421L, 421L, 
                      421L, 421L, 423L, 423L, 423L, 424L, 424L, 424L, 424L, 424L, 424L, 
                      506L, 506L, 506L, 562L, 562L, 562L, 562L, 562L, 562L, 562L, 562L, 
                      562L, 562L, 562L, 562L, 562L, 649L, 649L, 747L, 747L, 747L, 747L, 
                      747L, 747L, 798L, 866L, 866L, 866L, 866L, 866L, 930L, 930L, 930L, 
                      930L, 930L, 930L, 930L, 930L, 930L, 930L, 930L, 963L, 963L, 963L, 
                      963L, 963L, 963L, 963L, 963L, 1016L, 1016L, 1016L, 1016L, 1016L, 
                      1016L, 1016L, 1016L, 1066L, 1066L, 1102L, 1102L, 1102L, 1102L, 
                      1102L, 1185L, 1185L, 1185L, 1185L, 1185L, 1185L, 1185L, 1185L, 
                      1185L, 1185L, 1185L, 1185L, 1186L, 1186L, 1186L, 1186L, 1186L, 
                      1186L, 1186L, 1194L, 1194L, 1194L, 1194L, 1435L, 1531L, 1531L, 
                      1531L, 1531L, 1531L, 1531L, 1531L, 1547L, 1559L, 1559L, 1559L, 
                      1559L, 1559L), 
           dry.weight = c(24L, 116L, 52L, 30L, 142L, 40L, 
                      34L, 10L, 52L, 26L, 44L, 48L, 10L, 56L, 40L, 38L, 46L, 36L, 14L, 
                      24L, 130L, 34L, 24L, 56L, 30L, 28L, 52L, 386L, 46L, 46L, 16L, 
                      28L, 32L, 28L, 22L, 28L, 22L, 58L, 14L, 40L, 14L, 96L, 142L, 
                      114L, 46L, 34L, 46L, 114L, 130L, 38L, 134L, 44L, 42L, 26L, 34L, 
                      42L, 18L, 10L, 40L, 102L, 56L, 24L, 12L, 44L, 46L, 18L, 30L, 
                      52L, 58L, 52L, 4L, 64L, 14L, 74L, 206L, 30L, 108L, 20L, 46L, 
                      6L, 40L, 46L, 28L, 32L, 102L, 68L, 58L, 48L, 32L, 74L, 32L, 114L, 
                      58L, 32L, 28L, 48L, 6L, 32L, 26L, 64L, 108L, 34L, 46L, 84L, 28L, 
                      84L, 34L, 88L, 20L, 46L, 66L, 152L, 164L, 48L, 84L, 470L, 70L, 
                      42L, 294L, 110L, 174L, 126L, 54L, 872L, 48L, 312L, 62L, 162L, 
                      44L, 46L, 90L, 34L, 228L, 188L, 78L, 406L, 170L, 168L, 36L, 36L, 
                      76L, 24L, 30L, 58L, 82L, 124L, 32L, 76L, 36L, 88L, 94L, 26L), 
           healthy.pod = c(1L, 2L, 1L, 1L, 5L, 2L, 1L, 1L, 1L, 1L, 2L, 
                           1L, 2L, 3L, 2L, 1L, 1L, 1L, 2L, 1L, 6L, 1L, 1L, 1L, 1L, 1L, 
                           2L, 17L, 2L, 5L, 2L, 1L, 1L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 5L, 14L, 8L, 2L, 1L, 2L, 6L, 5L, 2L, 8L, 2L, 1L, 1L, 
                           2L, 1L, 1L, 1L, 1L, 4L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 2L, 
                           1L, 1L, 1L, 1L, 2L, 7L, 1L, 3L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 
                           2L, 2L, 1L, 1L, 1L, 2L, 1L, 3L, 4L, 1L, 2L, 2L, 1L, 3L, 1L, 
                           2L, 8L, 1L, 2L, 4L, 1L, 2L, 2L, 3L, 1L, 1L, 1L, 3L, 5L, 2L, 
                           3L, 12L, 2L, 2L, 7L, 2L, 4L, 3L, 1L, 19L, 1L, 6L, 1L, 4L, 
                           1L, 1L, 2L, 1L, 5L, 3L, 2L, 13L, 3L, 4L, 1L, 1L, 2L, 2L, 
                           1L, 4L, 2L, 3L, 1L, 1L, 1L, 3L, 2L, 2L), 
           infected.pods = c(0L, 
                           0L, 0L, 4L, 2L, 2L, 0L, 0L, 0L, 0L, 2L, 0L, NA, 1L, 0L, 0L, 
                           0L, 0L, 0L, NA, 0L, 0L, NA, 0L, 0L, 0L, 2L, 0L, 0L, 1L, 0L, 
                           0L, 0L, 3L, 0L, 2L, 0L, 2L, NA, NA, 0L, 0L, 2L, 0L, 0L, 0L, 
                           2L, 1L, 2L, 3L, 4L, 0L, 2L, NA, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
                           0L, 0L, NA, 0L, 0L, NA, 0L, 0L, NA, NA, 0L, 0L, 0L, 1L, 0L, 
                           0L, 0L, 0L, 0L, NA, NA, 2L, 2L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 
                           0L, NA, 2L, 0L, 2L, 0L, 1L, 0L, NA, 0L, 0L, 0L, 1L, NA, NA, 
                           0L, 0L, 0L, 0L, 0L, NA, 0L, 1L, 2L, 6L, 0L, 0L, 0L, 0L, NA, 
                           1L, 0L, 0L, 0L, 0L, 1L, 1L, NA, 0L, 0L, 0L, 0L, 0L, 0L, NA, 
                           0L, 0L, 0L, 0L, 0L, 0L, NA, 0L, 0L, 0L, NA, 0L, 0L, 0L, 2L, 
                           NA),
           date = structure(c(29L, 35L, 37L, 5L, 25L, 9L, 16L, 
                           13L, 33L, 7L, 11L, 11L, 8L, 3L, 19L, 7L, 12L, 17L, 4L, 36L, 
                           16L, 33L, 36L, 23L, 35L, 12L, 5L, 16L, 13L, 33L, 31L, 24L, 
                           9L, 37L, 17L, 9L, 35L, 11L, 36L, 8L, 33L, 7L, 25L, 29L, 28L, 
                           28L, 16L, 36L, 5L, 13L, 9L, 27L, 9L, 8L, 25L, 1L, 28L, 11L, 
                           35L, 16L, 22L, 5L, 29L, 36L, 34L, 31L, 8L, 5L, 17L, 36L, 
                           38L, 16L, 3L, 13L, 9L, 24L, 37L, 9L, 28L, 35L, 36L, 8L, 5L, 
                           23L, 28L, 9L, 34L, 37L, 33L, 12L, 15L, 35L, 36L, 15L, 17L, 
                           16L, 9L, 8L, 5L, 38L, 27L, 13L, 35L, 23L, 8L, 36L, 17L, 12L, 
                          17L, 24L, 24L, 36L, 35L, 5L, 16L, 16L, 33L, 25L, 29L, 28L, 
                           8L, 17L, 13L, 24L, 5L, 27L, 9L, 25L, 8L, 37L, 13L, 35L, 29L, 
                           28L, 25L, 36L, 28L, 29L, 30L, 30L, 17L, 26L, 39L, 37L, 24L, 
                           32L, 39L, 20L, 26L, 30L, 21L, 39L),
          .Label = c("02/09/2015", "03/08/2015", "04/07/2016", "04/08/2016", "04/08/2017", "04/09/2016", 
                           "05/05/2016", "05/10/2017", "06/07/2017", "06/10/2017", "07/04/2016", 
                          "07/04/2017", "07/06/2017", "07/07/2015", "07/09/2016", "07/09/2017", 
                          "07/10/2016", "08/01/2018", "08/06/2016", "08/06/2017", "08/08/2017", 
                          "08/10/2015", "09/05/2017", "09/12/2016", "10/03/2016", "10/05/2017", 
                          "10/11/2016", "11/01/2016", "11/02/2016", "11/09/2017", "11/11/2015", 
                          "11/11/2016", "12/01/2017", "14/12/2015", "16/03/2017", "16/11/2017", 
                          "17/02/2017", "18/12/2017", "20/11/2017"), class = "factor")), .Names = c("plotnb", 
                           "dry.weight", "healthy.pod", "infected.pods", "date"), row.names = c(286L, 
                           287L, 288L, 289L, 290L, 291L, 292L, 293L, 294L, 295L, 296L, 503L, 
                           924L, 925L, 926L, 927L, 928L, 929L, 930L, 931L, 932L, 933L, 934L, 
                            1365L, 1366L, 1367L, 1790L, 1791L, 1792L, 1846L, 1847L, 1848L, 
                             1849L, 1850L, 1981L, 1982L, 2366L, 2450L, 2451L, 2452L, 2453L, 
                             2454L, 2458L, 2459L, 2460L, 2461L, 2462L, 2463L, 2464L, 2465L, 
                            2466L, 2962L, 2963L, 2964L, 3212L, 3213L, 3214L, 3215L, 3216L, 
                           3217L, 3218L, 3219L, 3220L, 3221L, 3222L, 3223L, 3224L, 3531L, 
                            3532L, 3971L, 3972L, 3973L, 3974L, 3975L, 3976L, 4166L, 4387L, 
                            4388L, 4389L, 4390L, 4391L, 4605L, 4606L, 4607L, 4608L, 4609L, 
                            4610L, 4611L, 4612L, 4613L, 4614L, 4615L, 4747L, 4748L, 4749L, 
                            4750L, 4751L, 4752L, 4753L, 4754L, 5030L, 5031L, 5032L, 5033L, 
                            5034L, 5035L, 5036L, 5037L, 5252L, 5253L, 5411L, 5412L, 5413L, 
                            5414L, 5415L, 5761L, 5762L, 5763L, 5764L, 5765L, 5766L, 5767L, 
                            5768L, 5769L, 5770L, 5771L, 5772L, 5773L, 5774L, 5775L, 5776L, 
                            5777L, 5778L, 5779L, 5794L, 5795L, 5796L, 5797L, 6620L, 6807L, 
                            6808L, 6809L, 6810L, 6811L, 6812L, 6813L, 6840L, 6854L, 6855L, 
                            6856L, 6857L, 6858L), class = "data.frame")
monilia <- structure(list(plotnb = structure(c(24L, 24L, 24L, 156L, 162L, 
                                162L, 179L, 218L, 219L, 219L, 237L, 237L, 237L, 332L, 332L, 385L, 
                                385L, 385L), .Label = c("1", "10", "100", "101", "102", "103", 
                                                        "106", "107", "1073", "1074", "1078", "1079", "108", "1082", 
                                                        "1086", "1088", "109", "1091", "1097", "1098", "1099", "11", 
                                                        "1101", "1102", "1104", "1105", "1106", "1107", "1108", "1109", 
                                                        "111", "1112", "1116", "1118", "112", "1124", "1127", "1128", 
                                                        "1129", "1133", "1134", "1136", "1137", "1139", "1142", "1145", 
                                                        "1146", "1148", "115", "1151", "1152", "1153", "1155", "1157", 
                                                        "116", "1173", "118", "1180", "119", "1201", "121", "1242", "1243", 
                                                        "1248", "1259", "126", "1260", "1280", "1281", "1290", "1299", 
                                                        "13", "1302", "1318", "1334", "1347", "1365", "1375", "14", "1403", 
                                                        "1407", "1408", "141", "1412", "1445", "1446", "1447", "1451", 
                                                        "1452", "1453", "1455", "1457", "1467", "1472", "1476", "1483", 
                                                        "1492", "1519", "1524", "1525", "160", "172", "179", "18", "182", 
                                                        "183", "192", "2", "203", "21", "22", "220", "23", "235", "240", 
                                                        "247", "257", "258", "259", "26", "260", "261", "262", "264", 
                                                        "27", "271", "273", "274", "275", "276", "277", "278", "279", 
                                                        "28", "281", "283", "288", "290", "291", "292", "293", "295", 
                                                        "297", "298", "301", "302", "303", "305", "306", "307", "309", 
                                                        "31", "310", "313", "314", "318", "320", "321", "323", "324", 
                                                        "325", "326", "33", "331", "332", "334", "335", "336", "337", 
                                                        "34", "340", "341", "342", "343", "344", "346", "347", "348", 
                                                        "349", "350", "352", "355", "356", "357", "359", "360", "363", 
                                                        "364", "367", "368", "369", "372", "373", "375", "376", "377", 
                                                        "379", "38", "380", "381", "383", "386", "389", "391", "394", 
                                                        "398", "399", "4", "40", "400", "406", "41", "410", "411", "413", 
                                                        "414", "419", "423", "424", "430", "431", "433", "436", "437", 
                                                        "44", "441", "442", "443", "447", "453", "457", "46", "461", 
                                                        "470", "479", "486", "49", "490", "491", "497", "5", "50", "503", 
                                                        "52", "521", "524", "530", "533", "536", "539", "542", "547", 
                                                        "551", "552", "554", "556", "558", "561", "564", "568", "57", 
                                                        "574", "577", "579", "58", "580", "581", "582", "587", "588", 
                                                        "589", "590", "593", "597", "598", "599", "602", "604", "606", 
                                                        "607", "609", "61", "611", "613", "617", "62", "63", "637", "65", 
                                                        "655", "657", "66", "661", "662", "664", "666", "668", "67", 
                                                        "671", "683", "684", "685", "686", "688", "696", "698", "7", 
                                                        "702", "703", "704", "706", "71", "710", "711", "712", "717", 
                                                        "718", "721", "722", "724", "726", "733", "734", "735", "737", 
                                                        "739", "74", "740", "742", "743", "746", "747", "748", "755", 
                                                        "758", "76", "761", "762", "77", "773", "774", "777", "778", 
                                                        "78", "781", "783", "786", "789", "796", "797", "80", "801", 
                                                        "803", "804", "807", "808", "815", "816", "817", "819", "82", 
                                                        "820", "822", "823", "824", "826", "827", "828", "83", "830", 
                                                        "831", "836", "837", "838", "84", "844", "852", "853", "856", 
                                                        "858", "86", "863", "864", "865", "866", "87", "872", "876", 
                                                        "877", "880", "92", "94", "98", "99", "RDT53"), class = "factor"), 
           Fecha.Calificacion = structure(c(3L, 3L, 3L, 6L, 5L, 10L, 
                                            5L, 9L, 8L, 8L, 14L, 14L, 14L, 16L, 16L, 12L, 12L, 12L), .Label = c("04/07/2017", 
                                                                                                                "04/10/2017", "12/09/2017", "13/11/2017", "15/08/2017", "16/05/2017", 
                                                                                                                "17/05/2017", "18/05/2017", "20/06/2017", "20/09/2017", "20/12/2017", 
                                                                                                                "23/08/2017", "23/10/2017", "25/07/2017", "26/06/2017", "28/06/2017", 
                                                                                                                "28/11/2017", "29/11/2017", "30/05/2017", "31/05/2017"), class = "factor")), .Names = c("plotnb", 
                                                                                                                                                                                                        "Fecha.Calificacion"), row.names = c(59L, 60L, 61L, 400L, 412L, 
                                                                                                                                                                                                                                             413L, 456L, 552L, 553L, 554L, 591L, 592L, 593L, 768L, 769L, 907L, 
                                                                                                                                                                                                                                             908L, 909L), class = "data.frame")
phytophtora <-structure(list(plotnb = structure(c(17L, 17L, 17L, 17L, 80L, 
                                80L, 80L), .Label = c("1072", "1073", "1074", "1075", "1078", 
                                                      "1082", "1086", "1087", "1088", "1091", "1093", "1097", "1098", 
                                                      "1099", "1100", "1101", "1102", "1104", "1106", "1108", "1109", 
                                                      "1112", "1116", "1122", "1127", "1128", "1129", "1130", "1131", 
                                                      "1136", "1138", "1139", "1141", "1142", "1143", "1144", "1146", 
                                                      "1148", "1150", "1151", "1153", "1154", "1157", "1159", "375", 
                                                      "777", "778", "779", "781", "783", "788", "796", "799", "801", 
                                                      "803", "804", "807", "809", "812", "816", "819", "820", "823", 
                                                      "824", "827", "828", "836", "837", "838", "842", "843", "845", 
                                                      "846", "856", "858", "859", "861", "863", "864", "866", "867", 
                                                      "869", "871", "872", "875", "877", "RDT10"), class = "factor"), 
           Fecha.Calificacion = structure(c(3L, 7L, 1L, 2L, 7L, 7L, 
                                            7L), .Label = c("08/05/2017", "10/04/2017", "14/08/2017", 
                                                            "15/09/2017", "16/11/2017", "25/01/2018", "29/06/2017"), class = "factor")), .Names = c("plotnb", 
                                                            "Fecha.Calificacion"), row.names = c(36L, 37L, 38L, 39L, 170L, 
                                                            171L, 172L), class = "data.frame")

CODE

## using a sys.time() to check time taken.
start.time <- Sys.time()
## creating a column year to better group data.
yield$year <- "2015-2016"
yield$year[as.POSIXct(as.character(yield$date),format="%d/%m/%Y") >= "2017-07-01"] <- "2017-2018"
yield$year[as.POSIXct(as.character(yield$date),format="%d/%m/%Y") >= "2016-07-01" &
         as.POSIXct(as.character(yield$date),format="%d/%m/%Y") < "2017-07-01"] <- "2016-2017"
yield$year <- as.factor(yield$year)

## grouping data and calculating Total pod/ Total Healthy for each year.
yield.year <- c()
for (parcelle in unique(yield$plotnb)){
  for (year in levels(yield$year)){
    yield.subset <- yield[yield$plotnb==parcelle & yield$year==year,]
    TotalSane= sum(yield.subset$healthy.pod)
    TotalInfected= sum(yield.subset$infected.pods)
    TotalPreleved=0
    if (year=="2016-2017"){
      monilia.subset= monilia[monilia$plotnb==parcelle & as.POSIXct(as.character(monilia$Fecha.inoculacion),format="%d/%m/%Y") < "2017-07-01"  ,]
      phytophtora.subset= phytophtora[phytophtora$plotnb==parcelle & as.POSIXct(as.character(phytophtora$Fecha.inoculacion),format="%d/%m/%Y") < "2017-07-01" ,]
      TotalPreleved = nrow(monilia.subset)+ nrow(phytophtora.subset)
    } else if (year=="2017-2018"){
      monilia.subset= monilia[monilia$plotnb==parcelle & as.POSIXct(as.character(monilia$Fecha.inoculacion),format="%d/%m/%Y") >= "2017-07-01"  ,]
      phytophtora.subset= phytophtora[phytophtora$plotnb==parcelle & as.POSIXct(as.character(phytophtora$Fecha.inoculacion),format="%d/%m/%Y") >= "2017-07-01" ,]
      TotalPreleved = nrow(monilia.subset)+ nrow(phytophtora.subset)
    }
    TotalPod= TotalSane + TotalInfected + TotalPreleved
    TotalWeight= sum(yield.subset$dry.weight)
    yield.year <- rbind(yield.year,c(parcelle,TotalSane,TotalInfected,TotalPreleved,TotalPod,TotalWeight,year))
  }
}

## formating 
yield.year <- as.data.frame(yield.year)
colnames(yield.year) <- c("plotnb","TotalSane","TotalInfected","TotalPreleved","TotalPod","TotalWeight","year")

## calculating other data
yield.year[,"potential_yield"] <- as.numeric(as.vector(yield.year$TotalPod))*as.numeric(as.vector(yield.year$TotalWeight))/as.numeric(as.vector(yield.year$TotalSane))
yield.year[,"weight_per_pod"] <- as.numeric(as.vector(yield.year$TotalWeight))/as.numeric(as.vector(yield.year$TotalSane))
yield.year[,"TotalWeight"] <- as.numeric(as.character(yield.year$TotalWeight))
yield.year[,"TotalPod"] <- as.numeric(as.character(yield.year$TotalPod))

 ## cleaning
rm(list=c("monilia.subset","phytophtora.subset","yield.subset",
      "parcelle","TotalInfected","TotalPod","TotalPreleved","TotalSane","TotalWeight","year"))

## getting time.taken
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken
\$\endgroup\$
4
  • \$\begingroup\$ This line of code can't be right TotalPreleved = nrow(monilia.subset)+ nrow(phytophtora.subset) as the second data set is specified later if at all. So I skipped this task completely ou in my answer. Please revise your question how to summarise monilia and phytophtora. Concrete how to sum up for plotnb == 49 the Fecha.Calificacion. By the way it is three times the same date. \$\endgroup\$
    – Roman
    Commented Jul 18, 2018 at 14:35
  • \$\begingroup\$ @jimbou thanks for your review. Did you have an error with that codeline? The two datasets are normally specified just before doing the calculation. I am not sure to understand the second part of your comment about revision I need to do. There could indeed exist three observation of the same tree on the same date (one observation is a fruit) \$\endgroup\$
    – Untitpoi
    Commented Jul 18, 2018 at 14:46
  • \$\begingroup\$ Sorry you are right. So you only want to sum up the rows per plotnb and Year and add this information. \$\endgroup\$
    – Roman
    Commented Jul 18, 2018 at 15:30
  • \$\begingroup\$ can you comment why the monilia.subsetis not 3 for year=="2017-2018" and plotnb==49? \$\endgroup\$
    – Roman
    Commented Jul 18, 2018 at 15:47

2 Answers 2

2
\$\begingroup\$

Your date variables are stored as factors and everytime you need to work with them you have to transform them into dates. This needs time and can be avoided by transforming them once in the beginning.
This also makes your code easier readable.

You are calculating the year groups for each data.frame. As this is repetitive code you can wrap it in a function. Also you cann pull out the calculations from the for loop.

You are probably reading your data from a csv-file with read.table or your are using as.data.frame at some point. Setting stringsAsFactors = FALSE will make your life easier with less factors at the wrong places.

Another point is that you are growing the yield.year variable in each iteration of the for-loops. You can speed this up by initializing the variable with the final length before the looping.

You store the results in your for-loop in a matrix and thus everything is coerce to the same datatype. Using a data.frame for storing will let you use different data types and you don't have to use as.numeric later.

There are some missing values NA in your data. By using sum you will get NA as a result, if one summand is NA. To prevent this you should use na.rm = T.

Putting these things together I change your code and already got speed things up remarkably (roughly half the time it takes.)

To better compare different version I use the microbenchmark-package and wrapped the different versions in functions. I named your version old_version. (I changed Fecha.inoculacion to Fecha.Calificacion to match your data.)

This is the first improvment:

pre_initialized <- function(yield, monilia, phytophtora) {

  year_groups <- function(date) {
    factor(ifelse(date >= "2017-07-01",
                  "2017-2018",
                  ifelse(date >= "2016-07-01" & date < "2017-07-01",
                         "2016-2017",
                         "2015-2016")),
           levels = c("2015-2016", "2016-2017", "2017-2018")
    )
  }

  ## Convert date columns in correct format
  yield$date <- as.POSIXct(as.character(yield$date), format = "%d/%m/%Y")
  monilia$Fecha.Calificacion <- as.POSIXct(as.character(monilia$Fecha.Calificacion), format = "%d/%m/%Y")
  phytophtora$Fecha.Calificacion <- as.POSIXct(as.character(phytophtora$Fecha.Calificacion), format = "%d/%m/%Y")

  ## Set year intervals
  yield$year <- year_groups(yield$date)
  monilia$year <- year_groups(monilia$Fecha.Calificacion)
  phytophtora$year <- year_groups(phytophtora$Fecha.Calificacion)


  ## Initialize data.frame
  n_rows <- length(unique(yield$plotnb)) * length(levels(yield$year))
  yield.year <-
    data.frame(plotnb = rep(unique(yield$plotnb), each = length(levels(yield$year))),
               year = rep(levels(yield$year), times = length(unique(yield$plotnb))),
               TotalSane = numeric(n_rows),
               TotalInfected = numeric(n_rows),
               TotalPreleved = numeric(n_rows),
               TotalPod = numeric(n_rows),
               TotalWeight = numeric(n_rows),
               stringsAsFactors = F)


  ## grouping data and calculating Total pod/ Total Healthy for each year.
  for (parcelle in unique(yield$plotnb)) {
    for (year in levels(yield$year)) {
      yield.subset <- yield[yield$plotnb == parcelle & yield$year == year,]
      TotalSane = sum(yield.subset$healthy.pod)
      TotalInfected = sum(yield.subset$infected.pods)
      TotalPreleved = 0
      if (!year %in% c("2015-2016")) {
        monilia.subset = monilia[monilia$plotnb == parcelle & monilia$year == year, ]
        phytophtora.subset = phytophtora[phytophtora$plotnb == parcelle & phytophtora$year == year, ]
        TotalPreleved = nrow(monilia.subset) + nrow(phytophtora.subset)
      }
      TotalPod = TotalSane + TotalInfected + TotalPreleved
      TotalWeight = sum(yield.subset$dry.weight)
      yield.year[yield.year$plotnb == parcelle & yield.year$year == year,
                 c("TotalSane", "TotalInfected", "TotalPreleved", "TotalPod", "TotalWeight")] <-
        c(TotalSane, TotalInfected, TotalPreleved, TotalPod, TotalWeight)
    }
  }

  ## calculating other data
  yield.year[, "potential_yield"] <-
    yield.year$TotalPod * yield.year$TotalWeight / yield.year$TotalSane
  yield.year[, "weight_per_pod"] <-
    yield.year$TotalWeight / yield.year$TotalSane

  return(yield.year)
}

The main calculation is done in two nested for-loops. As you basically only group by plotnb and year this can be simplified by using the dplyr-package. Further I use tidyr for the functions complete (optain all combinations of plotnb and year) and replace_na.

This speeds up the calculations again.

library(dplyr)
library(tidyr)

dplyr_sol <- function(yield, monilia, phytophtora){

  year_groups <- function(date) {
    factor(case_when(date >= "2017-07-01" ~ "2017-2018", 
                     date >= "2016-07-01" & date < "2017-07-01" ~ "2016-2017", 
                     TRUE ~ "2015-2016"), 
           levels = c("2015-2016", "2016-2017", "2017-2018"))
  }

  yield <- yield %>% 
    mutate(date = as.Date(date, format = "%d/%m/%Y"), 
           year = year_groups(date)) %>% 
    complete(plotnb, year, 
             fill = list(dry.weight = 0, healthy.pod = 0, infected.pods = 0))

  monilia_sum <- monilia %>% 
    mutate(plotnb = as.integer(as.character(plotnb)),
           Fecha.Calificacion = as.Date(Fecha.Calificacion, format = "%d/%m/%Y"), 
           year = year_groups(Fecha.Calificacion)) %>% 
    group_by(plotnb, year) %>% 
    summarise(MonPreleved = n()) 

  phytophtora_sum <- phytophtora %>% 
    mutate(plotnb = as.integer(as.character(plotnb)),
           Fecha.Calificacion = as.Date(Fecha.Calificacion, format = "%d/%m/%Y"), 
           year = year_groups(Fecha.Calificacion)) %>% 
    group_by(plotnb, year) %>% 
    summarise(PhyPreleved = n()) 

  yield_sum <- yield %>% 
    group_by(plotnb, year) %>% 
    summarise(TotalSane = sum(healthy.pod),
              TotalInfected = sum(infected.pods),
              TotalWeight = sum(dry.weight))

  yield_year <- yield_sum %>% 
    left_join(monilia_sum, by = c("plotnb", "year")) %>% 
    left_join(phytophtora_sum, by = c("plotnb", "year")) %>% 
    replace_na(list(MonPreleved = 0, PhyPreleved = 0)) %>%  
    mutate(TotalPreleved = ifelse(!year == "2015-2016", MonPreleved + PhyPreleved, 0),
           TotalPod = TotalSane + TotalInfected + TotalPreleved, 
           potential_yield = TotalPod*TotalWeight/TotalSane, 
           weight_per_pod = TotalWeight/TotalSane) %>% 
    select(plotnb, TotalSane, TotalInfected, TotalPreleved, TotalPod, TotalWeight, 
           year, potential_yield, weight_per_pod)

  return(yield_year)
}

Here is the comparison of the different solutions:

library(microbenchmark)

microbenchmark(old_version(yield, monilia, phytophtora), 
               pre_initialized(yield, monilia, phytophtora), 
               dplyr_sol(yield, monilia, phytophtora))

Unit: milliseconds
                                         expr      min        lq      mean    median        uq       max neval cld
     old_version(yield, monilia, phytophtora) 98.24141 101.14597 105.26297 102.84971 105.31294 207.91746   100   c
 pre_initialized(yield, monilia, phytophtora) 47.40690  49.17982  52.00089  50.65794  52.80491 117.60703   100  b 
       dplyr_sol(yield, monilia, phytophtora) 19.63437  20.70025  22.49237  21.29610  22.32720  98.36952   100 a  
\$\endgroup\$
2
\$\begingroup\$

For the first part you can try a tidyverse

library(tidyverse)
yield %>% 
  as.tibble() %>%
  mutate(date=as.POSIXct(date,format="%d/%m/%Y")) %>%
  mutate(Year=cut(date, as.POSIXct(seq.Date(as.Date("2015-07-01"), by="year", length.out = 4)))) %>% 
  mutate(Year=factor(Year, labels = c("2015-2016","2016-2017","2017-2018"))) %>% 
  group_by(plotnb, Year) %>% 
  summarise_at(vars(dry.weight, healthy.pod, infected.pods),
               funs(Total=sum(., na.rm=T))) %>% 
  complete(plotnb, Year, fill = list(dry.weight_Total = 0, 
                                     healthy.pod_Total = 0,infected.pods_Total=0)) %>% 
  mutate(TotalPreleved=0) %>% 
  mutate(TotalPod=healthy.pod_Total + infected.pods_Total + TotalPreleved)
# A tibble: 90 x 7
# Groups:   plotnb [30]
   plotnb Year      dry.weight_Total healthy.pod_Total infected.pods_Total TotalPreleved TotalPod
    <int> <fct>                <dbl>             <dbl>               <dbl>         <dbl>    <dbl>
 1     49 2015-2016              236                 9                   4             0       13
 2     49 2016-2017              230                 5                   0             0        5
 3     49 2017-2018              104                 4                   6             0       10
 4     89 2015-2016               48                 1                   0             0        1
 5     89 2016-2017                0                 0                   0             0        0
 6     89 2017-2018                0                 0                   0             0        0
 7    158 2015-2016               78                 3                   0             0        3
 8    158 2016-2017              152                 7                   1             0        8
 9    158 2017-2018              164                 9                   0             0        9
10    159 2015-2016                0                 0                   0             0        0
# ... with 80 more rows
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.