Example 1: photosystem efficiency in Arabidopsis

The used in this vignette contains data from an experiment in the Phenovator platform (WUR, Netherlands, (Flood et al. 2016)) with Arabidopsis plants. It consists of one experiment with 1440 plants grown in a growth chamber with different light intensity. The data set called “PhenovatorDat1” is included in the package.

The number of tested genotypes (Genotype) is 192 with 6-7 replicates per genotype (Replicate). Four reference genotypes were also tested with 15 or 30 replicates. The studied trait is the photosystem II efficiency (EffpsII) extracted from the pictures over time (Rooijen et al. 2017). The unique ID of the plant is recorded (pos), together with the pot position in row (x) and in column (y). The data set also includes factors from the design: the position of the camera (Image_pos) and the pots table (Basin).

data("PhenovatorDat1")

Layout plot

The first type of plot displays the layout of the experiment as a grid using the row and column coordinates. The default option creates plots of all time points in the TP object. This can be restricted to a selection of time points using their number in the option timePoints. If repId was specified when creating the TP object, replicate blocks are delineated by a black line. Missing plots are indicated in white enclosed with a bold black line. This type of plot allows checking the design of the experiment.

## Plot the layout for the third time point.
plot(phenoTP, 
     plotType = "layout",
     timePoints = 3)

Here, the third time point is displayed which corresponds to the 1^st of June 2018 at 11:37. Note that the title can be manually changed using the title option. This plot can be extended by highlighting interesting genotypes in the layout. Hereafter the check genotypes are highlighted:

## Plot the layout for the third time point with the check genotypes highlighted.
plot(phenoTP, 
     plotType = "layout",
     timePoints = 3,  
     highlight = c("check1", "check2", "check3", "check4"))

It is possible to add the labels of the genotypes to the layout.

## Plot the layout for the third time point.
plot(phenoTP, 
     plotType = "layout",
     timePoints = 3,  
     highlight = c("check1", "check2", "check3", "check4"),
     showGeno = TRUE)

We can visualize the raw data of a given trait on the layout, as a heatmap. This type of plot gives a first indication of the spatial variability at a given time point. This can be further investigated with the spatial modeling.

## Plot the layout for the third time point.
plot(phenoTP, 
     plotType = "layout",
     timePoints = 3,  
     traits = "EffpsII")

Raw data plot

Raw data can be displayed per genotype with one color per plotId.

## Create the raw data time courses for three genotypes.
plot(phenoTP, 
     traits = "EffpsII",
     plotType = "raw",
     genotypes = c("G001", "G002", "check1"))

Boxplot

Boxplots can be made to visually assess the variability of the trait(s) in the TP object. By default a box is plotted per time point for the specified trait using all time points.

## Create a boxplot for "EffpsII" using the default all time points.
plot(phenoTP, 
     plotType = "box",
     traits = "EffpsII") 

Correlation plot

Finally, a plot of the correlations between the observations in time for a specified trait can be made. The order of the plot is chronological and by default all time points are used.

## Create a correlation plot for "EffpsII" for a selection of time points.
plot(phenoTP, 
     plotType = "cor",
     traits = "EffpsII",
     timePoints = seq(from = 1, to = 73, by = 5))

Spatial model using SpATS

When SpATS is used for modeling, an extra spatial term is included in the model. This spatial component is composed using the PSANOVA function in the SpATS package which uses 2-dimensional smoothing with P-splines as described in Lee, Durbán, and Eilers (2013) and in Rodríguez-Álvarez et al. (2018).

In the example below the model is fitted for a few points.

## Fit a model for a few time points.
modPhenoSp <- fitModels(TP = phenoTPOut, 
                        trait = "EffpsII",
                        timePoints = c(1, 33, 36, 54, 73))
summary(modPhenoSp)
#> Models in modPhenoSp where fitted for experiment Phenovator.
#> 
#> It contains 5 time points.
#> The models were fitted using SpATS.

Model plots

Several types of plots can be made for fitted models. A spatial plot consists of spatial plots of the raw data, fitted values, residuals and either BLUEs or BLUPs, and a histogram of the BLUEs or BLUPs. When SpATS is used for modeling an extra plot with the fitted spatial trend is included.

plot(modPhenoSp,
     timePoints = 36,
     plotType = "spatial",
     spaTrend = "percentage")

It is also possible to create a time lapse of the ratio of spatial trends over time. The scale is the same as previously described. The time lapse is always written to an output file.

plot(modPhenoSp, 
     plotType = "timeLapse",
     outFile = "TimeLapse_modPhenoSp.gif")

In the rawPred plot the raw data (colored dots, one color per plotId) is plotted overlaid with the predicted values (black dots) from the fitted model.

Note: The model has been fitted on a limited number of time points, which is reflected in the plots. For plots using more time point see this vignette.

plot(modPhenoSp, 
     plotType = "rawPred",
     genotypes = c("G007", "G058")) 

In the corrPred plot the corrected data (colored dots, one color per plotId) is plotted overlaid with the predicted values (black dots) from the fitted model.

plot(modPhenoSp, 
     plotType = "corrPred",
     genotypes = c("G007", "G058"))

The last three types of plot display different model parameters over time. Plot type herit plots the heritability over time.

plot(modPhenoSp, 
     plotType = "herit",
     yLim = c(0.5, 1))

Plot type variance plots the residual, column and row variances over time. These plots can serve as diagnostics of the experiment.

plot(modPhenoSp, 
     plotType = "variance")

Plot type effDim plots the effective dimension from models fitted using SpATS over time.

plot(modPhenoSp, 
     plotType = "effDim",
     whichED = c("colId", "rowId", "fColRow","colfRow", "surface"),
     EDType = "ratio")

Extracting model results

All results that can be extracted are shown in the table below. The first column contains the function names. The second column contains a short description of the result that will be extracted and, where needed, states for which modeling engines it can be extracted.

## Extract the genotypic predictions for one time point: 
genoPredSp <- getGenoPred(modPhenoSp)

Model fitting

To fit the psHDM model, we use the fitSplineHDM() function with cubic (bdeg \(= 3\)) B-spline basis functions with nseg \(= 20\) segments, and second order penalties (pord\(= 2\)) to represent functions at the three levels of the hierarchy (population, genotype and plot).

## Fit P-spline HDM.
fit.psHDM  <- fitSplineHDM(inDat = spatCorrectedVator,
                           genotypes = subGenoVator,
                           trait = "EffpsII_corr",
                           useTimeNumber = TRUE,
                           timeNumber = "timeNumHour",                           
                           pop = "pop",
                           genotype = "genotype",
                           plotId = "plotId",
                           weights = "wt",
                           difVar = list(geno = FALSE, plot = FALSE),
                           smoothPop = list(nseg = 20, bdeg = 3, pord = 2),
                           smoothGeno = list(nseg = 20, bdeg = 3, pord = 2),
                           smoothPlot = list(nseg = 20, bdeg = 3, pord = 2),
                           trace = FALSE)

The resulting object, in this case fit.psHDM, contains different information about the data structure, the fitting process, and three data frames with the estimated curves at each of the three-levels of the hierarchy (explore fit.psHDM$popLevel, fit.psHDM$genoLevel, and fit.psHDM$plotLevel). That is, it contains the estimated trajectories and deviations, as well as their first and second-order derivatives. For a detailed description of the returned values see help(fitSplineHDM). All these curves can be used as input to extract time‐independent parameters to characterise genotypes (as illustrated in the final section of this vignette).

Model prediction

The predict.psHDM() function can be used to obtain predictions from a fitted psHDM model (obtained using the fitSplineHDM() function; see above). In particular, this function allows obtaining predictions (estimated curves at each level of the hierarchy) on a dense grid of time points. Also, it allows the calculation of standard errors. These standard errors can be used to construct (approximate) pointwise confidence intervals for the estimated curves.

In the following code, we use the fit.psHDM object to make predictions at the three levels of the hierarchy (pred = list(pop = TRUE, geno = TRUE, plot = TRUE)), and to obtain standard errors at the population and genotype levels (se = list(pop = TRUE, geno = TRUE, plot = FALSE)). The original data is measured at \(73\) time points, but predictions are obtained at \(100\) time points in the same range than the original time points (argument newtimes).

## Predict P-spline HDM.
pred.psHDM <- predict(object = fit.psHDM,
                      newtimes = seq(min(fit.psHDM$time[["timeNumber"]]),
                                     max(fit.psHDM$time[["timeNumber"]]),
                                     length.out = 100),
                      pred = list(pop = TRUE, geno = TRUE, plot = TRUE),
                      se = list(pop = TRUE, geno = TRUE, plot = FALSE),
                      trace = FALSE)

As result, three data frames with predictions (and standard errors) at population (pred.psHDM$popLevel), genotype (pred.psHDM$GenoLevel) and plot (pred.psHDM$plotLevel) levels are returned.

Model plot

The plot.psHDM() function plots psHDM objects. We note that objects of class psHDM can be obtained using both fitSplineHDM() and predict.psHDM() functions. In both cases, the resulting object contains information about estimated trajectories, deviations and first-order derivatives at the three levels of the hierarchy. As such, plots of these curves can be obtained. In addition, when plots are obtained from an object obtained using the predict.psHDM() function, \(95\%\) pointwise confidence intervals are also depicted.

To illustrate the usage of function plot.psHDM(), we use here the object pred.psHDM obtained in the prediction section.

• If plotType = "popTra", the estimated population trajectory (\(\hat{f}_p(t)\)) is depicted and its \(95\%\) pointwise confidence intervals. Additionally, the grey lines represent the observed trait that is used in the fitSplineHDM function.

plot(pred.psHDM, plotType = "popTra", themeSizeHDM = 10)

• If plotType = "popGenoTra", estimated population (\(\hat{f}_p(t)\)) and genotype-especific (\(\hat{f}_p(t)+\hat{f}_{pg}(t)\)) trajectories are depicted for all genotypes. \(95\%\) pointwise confidence intervals are depicted for the population trajectory. We let the user explore the plot for the first-order derivative of the estimated population and genotype trajectories (plotType = "popGenoDeriv").

plot(pred.psHDM, plotType = "popGenoTra", themeSizeHDM = 10)

• If plotType = "GenoDev", estimated genotype-specific deviations (\(\hat{f}_{pg}(t)\)) are depicted for all genotypes.

plot(pred.psHDM, plotType = "genoDev", themeSizeHDM = 10)

Finally, if we are interested in obtaining plots at the plot level, we use plotType = "genoPlotTra". Here, estimated genotype (\(\hat{f}_p(t)+\hat{f}_{pg}(t)\)) and plot-specific (\(\hat{f}_p(t)+\hat{f}_{pg}(t)+\hat{f}_{pgi}(t)\)) trajectories are depicted for all plots separately for a selection of genotypes. Also, \(95\%\) pointwise confidence intervals are depicted for the estimated genotype-specific trajectories.

plot(pred.psHDM, 
     plotType = "genoPlotTra", 
     themeSizeHDM = 10)