tests/testthat/testRice.r

In imageData: Aids in Processing and Plotting Data from a Lemna-Tec Scananalyzer

#devtools::test("asremlPlus")
context("model_selection")

cat("#### Test using Rice baby example\n")
test_that("Rice2015_imageData", {
  skip_if_not_installed("imageData")
  skip_on_cran()
  library(dae)
  library(ggplot2)
  library(imageData)
  
  ## A dummy example to illustrate the use of imageData
  #'# Step 1: Import the raw data
  data(RiceRaw.dat)
  raw.dat <- RiceRaw.dat[1:280, ]
  raw.dat$Smarthouse <- 1
  
  #'# Step 2: Select imaging variables and add covariates and factors (produces longi.dat)
  longi.prime.dat <- longitudinalPrime(data=raw.dat, smarthouse.lev=1, 
                                       idcolumns = c("Genotype.ID", "Treatment.1", "Treatment.2"))
  
  longi.dat <- designFactors(longi.prime.dat, insertName = "xDays",
                             nzones = 1, nlanesperzone = 1, nmainplotsperlane = 10, 
                             designfactorMethod="StandardOrder")
  testthat::expect_equal(nrow(longi.dat), 280)
  testthat::expect_equal(ncol(longi.dat), 44)
  
  
  #'## Plot the imaging times for 20 carts from a single Lane
  testthat::expect_silent(imagetimesPlot(longi.dat, intervals = "Days", 
                                         timePositions = "Hour",
                                         ggplotFuncs=list(ggplot2::geom_line(aes(group=Snapshot.ID.Tag, 
                                                                                 colour=Lane)))))
  
  #'## Particular edits to longi.dat
  longi.dat <- within(longi.dat, 
                      { 
                        Days.after.Salting <- as.numfac(Days) - 29
                      })
  
  #'# Step 3: Form derived traits that result in a value for each observation
  #'### Set responses
  responses.image <- c("Area")
  responses.smooth <- paste(responses.image, "smooth", sep=".")
  
  #'## Form growth rates for each observation of a subset of responses by differencing
  longi.dat <- splitContGRdiff(longi.dat, responses.image, 
                               INDICES="Snapshot.ID.Tag",
                               which.rates = c("AGR","RGR"))
  testthat::expect_equal(nrow(longi.dat), 280)
  testthat::expect_equal(ncol(longi.dat), 48)
  
  #'## Form Area.WUE 
  longi.dat <- within(longi.dat, 
                      { 
                        Area.WUE <- WUI(Area.AGR*Days.diffs, Water.Loss)
                      })
  
  #'## Add cumulative responses 
  longi.dat <- within(longi.dat, 
                      { 
                        Water.Loss.Cum <- unlist(by(Water.Loss, Snapshot.ID.Tag, 
                                                    cumulate, exclude.1st=TRUE))
                        WUE.cum <- Area / Water.Loss.Cum 
                      })
  testthat::expect_equal(nrow(longi.dat), 280)
  testthat::expect_equal(ncol(longi.dat), 51)
  
  #'# Step 4: Fit splines to smooth the longitudinal trends in the primary traits and calculate their growth rates
  #'
  #'## Smooth responses
  #+
  for (response in c(responses.image, "Water.Loss"))
    longi.dat <- splitSplines(longi.dat, response, x="xDays", INDICES = "Snapshot.ID.Tag", 
                              df = 4)
  longi.dat <- with(longi.dat, longi.dat[order(Snapshot.ID.Tag, xDays), ])
  testthat::expect_equal(nrow(longi.dat), 280)
  testthat::expect_equal(ncol(longi.dat), 53)
  
  #'## Loop over smoothed responses, forming growth rates by differences
  #+
  responses.GR <- paste(responses.smooth, "AGR", sep=".")
  longi.dat <- splitContGRdiff(longi.dat, responses.smooth, 
                               INDICES="Snapshot.ID.Tag",
                               which.rates = c("AGR","RGR"))
  testthat::expect_equal(nrow(longi.dat), 280)
  testthat::expect_equal(ncol(longi.dat), 55)
  
  #'## Finalize longi.dat
  longi.dat <- with(longi.dat, longi.dat[order(Snapshot.ID.Tag, xDays), ])
  
  #'# Step 5: Do exploratory plots on unsmoothed and smoothed longitudinal data
  responses.longi <- c("Area","Area.AGR","Area.RGR", "Area.WUE")
  responses.smooth.plot <- c("Area.smooth","Area.smooth.AGR","Area.smooth.RGR")
  titles <- c("Total area (1000 pixels)", 
              "Total area AGR (1000 pixels per day)", "Total area RGR (per day)",
              "Total area WUE (1000 pixels per mL)")
  titles.smooth<-titles
  nresp <- length(responses.longi)
  limits <- list(c(0,1000), c(-50,125), c(-0.05,0.40), c(0,30))
  
  #' ### Plot unsmoothed profiles for all longitudinal  responses 
  #+ "01-ProfilesAll"
  klimit <- 0
  for (k in 1:nresp)
  { 
    klimit <- klimit + 1
    plt <- longiPlot(data = longi.dat, response = responses.longi[k], 
                     y.title = titles[k], x="xDays+35.42857143", printPlot=FALSE)
    plt <- plt + geom_vline(xintercept=29, linetype="longdash", size=1) +
      scale_x_continuous(breaks=seq(28, 42, by=2)) + 
      scale_y_continuous(limits=limits[[klimit]])
    if (k < 4)
      testthat::expect_silent(print(plt))
    else
      testthat::expect_warning(print(plt))
  }
  
  
  #' ### Plot smoothed profiles for all longitudinal responses - GRs by difference
  #+ "01-SmoothedProfilesAll"
  nresp.smooth <- length(responses.smooth.plot)
  limits <- list(c(0,1000), c(0,100), c(0.0,0.40))
  for (k in 1:nresp.smooth)
  { 
    plt <- longiPlot(data = longi.dat, response = responses.smooth.plot[k], 
                     y.title = titles.smooth[k], x="xDays+35.42857143", printPlot=FALSE)
    plt <- plt + geom_vline(xintercept=29, linetype="longdash", size=1) +
      scale_x_continuous(breaks=seq(28, 42, by=2)) + 
      scale_y_continuous(limits=limits[[k]])
    testthat::expect_silent(print(plt))
  }
  
  
  #'### AGR anomalies - plot without anomalous plants followed by plot of anomalous plants
  #+ "01-0254-AGRanomalies"
  anom.ID <- vector(mode = "character", length = 0L)
  response <- "Area.smooth.AGR"
  cols.output <- c("Snapshot.ID.Tag", "Smarthouse", "Lane", "Position", 
                   "Treatment.1", "Genotype.ID", "Days")
  anomalous <- anomPlot(longi.dat, response=response, lower=2.5, start.time=40, 
                        x = "xDays+35.42857143", vertical.line=29, breaks=seq(28, 42, by=2), 
                        whichPrint=c("innerPlot"), y.title=response)
  subs <- subset(anomalous$data, Area.smooth.AGR.anom & Days==42)
  if (nrow(subs) == 0)
  { cat("\n#### No anomalous data here\n\n")
  } else
  { 
    subs <- subs[order(subs["Smarthouse"],subs["Treatment.1"], subs[response]),]
    print(subs[c(cols.output, response)])
    anom.ID <- unique(c(anom.ID, subs$Snapshot.ID.Tag))
    outerPlot <- anomalous$outerPlot  + geom_text(data=subs,
                                                  aes_string(x = "xDays+35.42857143", 
                                                             y = response, 
                                                             label="Snapshot.ID.Tag"), 
                                                  size=3, hjust=0.7, vjust=0.5)
    print(outerPlot)
  }
  testthat::expect_equal(nrow(subs), 0)  
  
  #'# Step 6: Form single-value plant responses in Snapshot.ID.Tag order.
  #'
  #'## 6a) Set up a data frame with factors only
  #+
  cart.dat <- longi.dat[longi.dat$Days == 31, 
                        c("Smarthouse","Lane","Position","Snapshot.ID.Tag",
                          "xPosn","xMainPosn",
                          "Zones","xZones","SHZones","ZLane","ZMainplots", "Subplots",
                          "Genotype.ID","Treatment.1")]
  cart.dat <- cart.dat[do.call(order, cart.dat), ]
  
  #'## 6b) Get responses based on first and last date.
  #'
  #'### Observation for first and last date
  cart.dat <- cbind(cart.dat, getDates(responses.image, data = longi.dat, 
                                       which.times = c(31), suffix = "first"))
  cart.dat <- cbind(cart.dat, getDates(responses.image, data = longi.dat, 
                                       which.times = c(42), suffix = "last"))
  cart.dat <- cbind(cart.dat, getDates(c("WUE.cum"), 
                                       data = longi.dat, 
                                       which.times = c(42), suffix = "last"))
  responses.smooth <- paste(responses.image, "smooth", sep=".")
  cart.dat <- cbind(cart.dat, getDates(responses.smooth, data = longi.dat, 
                                       which.times = c(31), suffix = "first"))
  cart.dat <- cbind(cart.dat, getDates(responses.smooth, data = longi.dat, 
                                       which.times = c(42), suffix = "last"))
  testthat::expect_equal(nrow(cart.dat), 20)
  testthat::expect_equal(ncol(cart.dat), 19)
  
  #'### Growth rates over whole period.
  #+
  tottime <- 42 - 31
  cart.dat <- within(cart.dat, 
                     { 
                       Area.AGR <- (Area.last - Area.first)/tottime
                       Area.RGR <- log(Area.last / Area.first)/tottime
                     })
  
  #'### Calculate water index over whole period
  cart.dat <- merge(cart.dat, 
                    intervalWUI("Area", water.use = "Water.Loss", 
                                start.times = c(31), 
                                end.times = c(42), 
                                suffix = NULL, 
                                data = longi.dat, include.total.water = TRUE),
                    by = c("Snapshot.ID.Tag"))
  names(cart.dat)[match(c("Area.WUI","Water.Loss.Total"),names(cart.dat))] <- c("Area.Overall.WUE", 
                                                                                "Water.Loss.Overall")
  cart.dat$Water.Loss.rate.Overall <- cart.dat$Water.Loss.Overall / (42 - 31)
  testthat::expect_equal(nrow(cart.dat), 20)
  testthat::expect_equal(ncol(cart.dat), 25)
  
  #'## 6c) Add growth rates and water indices for intervals
  #'### Set up intervals
  #+
  start.days <- list(31,35,31,38)
  end.days <- list(35,38,38,42)
  suffices <- list("31to35","35to38","31to38","38to42")
  
  #'### Rates for specific intervals from the smoothed data by differencing
  #+
  for (r in responses.smooth)
  { for (k in 1:length(suffices))
  { 
    cart.dat <- merge(cart.dat, 
                      intervalGRdiff(r, 
                                     which.rates = c("AGR","RGR"), 
                                     start.times = start.days[k][[1]], 
                                     end.times = end.days[k][[1]], 
                                     suffix.interval = suffices[k][[1]], 
                                     data = longi.dat),
                      by = "Snapshot.ID.Tag")
  }
  }
  
  #'### Water indices for specific intervals from the unsmoothed and smoothed data
  #+
  for (k in 1:length(suffices))
  { 
    cart.dat <- merge(cart.dat, 
                      intervalWUI("Area", water.use = "Water.Loss", 
                                  start.times = start.days[k][[1]], 
                                  end.times = end.days[k][[1]], 
                                  suffix = suffices[k][[1]], 
                                  data = longi.dat, include.total.water = TRUE),
                      by = "Snapshot.ID.Tag")
    names(cart.dat)[match(paste("Area.WUI", suffices[k][[1]], sep="."), 
                          names(cart.dat))] <- paste("Area.WUE", suffices[k][[1]], sep=".")
    cart.dat[paste("Water.Loss.rate", suffices[k][[1]], sep=".")] <- 
      cart.dat[[paste("Water.Loss.Total", suffices[k][[1]], sep=".")]] / 
      ( end.days[k][[1]] - start.days[k][[1]])
  }
  
  cart.dat <- with(cart.dat, cart.dat[order(Snapshot.ID.Tag), ])
  testthat::expect_equal(nrow(cart.dat), 20)
  testthat::expect_equal(ncol(cart.dat), 49)
  
  #'# Step 7: Form continuous and interval SIITs
  #'
  #'## 7a) Calculate continuous
  #+
  cols.retained <-  c("Snapshot.ID.Tag","Smarthouse","Lane","Position",
                      "Days","Snapshot.Time.Stamp", "Hour", "xDays",
                      "Zones","xZones","SHZones","ZLane","ZMainplots",
                      "xMainPosn", "Genotype.ID")
  responses.GR <- c("Area.smooth.AGR","Area.smooth.AGR","Area.smooth.RGR")
  suffices.results <- c("diff", "SIIT", "SIIT")
  responses.SIIT <- unlist(Map(paste, responses.GR, suffices.results,sep="."))
  
  longi.SIIT.dat <- 
    twoLevelOpcreate(responses.GR, longi.dat, suffices.treatment=c("C","S"),
                     operations = c("-", "/", "/"), suffices.results = suffices.results, 
                     columns.retained = cols.retained, 
                     by = c("Smarthouse","Zones","ZMainplots","Days"))
  longi.SIIT.dat <- with(longi.SIIT.dat, 
                         longi.SIIT.dat[order(Smarthouse,Zones,ZMainplots,Days),])
  testthat::expect_equal(nrow(longi.SIIT.dat), 140)
  testthat::expect_equal(ncol(longi.SIIT.dat), 22)
  
  #' ### Plot SIIT profiles 
  #' 
  #+ "03-SIITProfiles"
  k <- 2
  nresp <- length(responses.SIIT)
  limits <- with(longi.SIIT.dat, list(c(min(Area.smooth.AGR.diff, na.rm=TRUE),
                                        max(Area.smooth.AGR.diff, na.rm=TRUE)),
                                      c(0,3),
                                      c(0,1.5)))
  #Plots
  for (k in 1:nresp)
  { 
    plt <- longiPlot(data = longi.SIIT.dat, x="xDays+35.42857143", 
                     response = responses.SIIT[k], 
                     y.title=responses.SIIT[k], 
                     facet.x="Smarthouse", facet.y=".", printPlot=FALSE, )
    plt <- plt + geom_vline(xintercept=29, linetype="longdash", size=1) +
      scale_x_continuous(breaks=seq(28, 42, by=2)) + 
      scale_y_continuous(limits=limits[[k]])
    testthat::expect_silent(print(plt))
  }
  
  #'## 7b) Calculate interval SIITs and check for large values for SIIT for Days 31to35
  #+ "01-SIITIntClean"
  suffices <- list("31to35","35to38","31to38","38to42")
  response <- "Area.smooth.RGR.31to35"
  SIIT <- paste(response, "SIIT", sep=".")
  responses.SIITinterval <- as.vector(outer("Area.smooth.RGR", suffices, paste, sep="."))
  
  cart.SIIT.dat <- twoLevelOpcreate(responses.SIITinterval, cart.dat,
                                    suffices.treatment=c("C","S"), 
                                    suffices.results="SIIT", 
                                    columns.suffixed="Snapshot.ID.Tag")
  testthat::expect_equal(nrow(cart.SIIT.dat), 10)
  testthat::expect_equal(ncol(cart.SIIT.dat), 23)
  tmp<-na.omit(cart.SIIT.dat)
  print(summary(tmp[SIIT]))
  big.SIIT <- with(tmp, tmp[tmp[SIIT] > 1.15, c("Snapshot.ID.Tag.C","Genotype.ID",
                                                paste(response,"C",sep="."), 
                                                paste(response,"S",sep="."), SIIT)])
  big.SIIT <- big.SIIT[order(big.SIIT[SIIT]),]
  testthat::expect_equal(nrow(big.SIIT), 0)
  testthat::expect_equal(ncol(big.SIIT), 5)
  print(big.SIIT)
  plt <- ggplot(tmp, aes_string(SIIT)) +
    geom_histogram(aes(y = ..density..), binwidth=0.05) +
    geom_vline(xintercept=1.15, linetype="longdash", size=1) +
    theme_bw() + facet_grid(Smarthouse ~.)
  print(plt)
  plt <- ggplot(tmp, aes_string(x="Smarthouse", y=SIIT)) +
    geom_boxplot() + theme_bw()
  testthat::expect_silent(print(plt))

})