@neo0351 we have scales listed as a depends. not sure I know where we use it?

apply_plsr(..) to take plsr coefficients/permutation coefficients and input spectra matrix to generate new values in a dataframe.

Options should be:
input coefficients as either R object or a specific file(s)
Input reflectance as either an R matrix or a specific file
wavelength prefix (e.g. "Wave_")
metadata_columns - those in the input spectra to keep in the output object

I don't see this overly complicated but we will want to have a check that wavelengths match. also need to indicate if there is an intercept and where it is.

Im sure other considerations will come up.

Use EcoSIS spectra and built-in coefficients as examples

As part of this build in the Serbin et al 2019 LMA coefficients

Using map/for.each like syntax to break the # of iterations into a set of separate parallel runs. Speed up for larger datasets

If ensuring the validation data set has a similar distribution as the calibration data set is important for model validation and is achieved by stratifying the random selection procedure used to split the data set, I think it should equally equally important for permuting the calibration data set to select the optimal number of components.
I suggest to improve find_optimal_components by allowing it to stratify the random sampling in each permutation.
In case you find it interesting, although I am far from en expert, I have already modified pls_permutation and find_optimal_components to allow for the input of a groups object of the form c("var1", "var2"..., "varn"). I names these functions as pls_permut_by_groups and find_optimal_comp_groups. pls_permut_by_groupsnow seems to perform well. I'm not sure about find_optimal_comp_groups, as I get a the following error:
Error: 'pls_permut_by_goups' is not an exported object from 'namespace:spectratrait'
Called from: getExportedValue(pkg, name)
I understand the issue, but I don't know how to tell getExportedValue that the new function is not in spectratrait.
Of course, this should not be an issue if you find the function could be included within spectratrait.

• Convert to package (shawn)
• Remove unused depends (@neo0351)
• Update README - revise links and update contains and other text (@neo0351 )
• Update example scripts in inst/scripts (shawn)
• Update vignettes (shawn)
• Add tests (all)
• TEST group testing sesh, squash bugs (all)
• Create internal dataset to use in an example not dependent on EcoSIS
• Create additional supplemental example that sources a package dataset (shawn)
• Update example script to include supplemental figure captions
  ...
• more here

Target var: Narea

> lapply(list.of.packages, library, character.only = TRUE)
[[1]]
 [1] "spectratrait" "gridExtra"    "ggplot2"      "plotrix"      "here"         "reshape2"    
 [7] "dplyr"        "pls"          "stats"        "graphics"     "grDevices"    "utils"       
[13] "datasets"     "methods"      "base"        

I don't think that's correct. Maybe at some point, I screwed up that syntax. clearly its wrong and it should be doing "library(package), ...)"

Responding to user feedback I noticed that some plots use index not wavelength on the X axis.
https://github.com/TESTgroup-BNL/spectratrait/blob/master/vignettes/ely_etal_ex1_files/figure-gfm/unnamed-chunk-12-1.png

This should be an easy fix

Aparently, the function performs a stratified random sampling for cal_data, but always takes the last rows of the dataset for val_data, even reusing cases already selected for cal_data:

plot<- rep(c("plot1", "plot2", "plot3"),each=42)
season<- rep(1:6, 21)
disease<- c(rep(0,84), rep(1,42))
d<- seq(1:126)
df <- data.frame(plot,season,disease,d)
df <- df %>% mutate(id=row_number())
names(df)
[1] "plot" "season" "disease" "d" "id"
df$id
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
[19] 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
[37] 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
[55] 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
[73] 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
[91] 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
[109] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

split_data <- spectratrait::create_data_split(dataset=df, approach="dplyr", split_seed=7529075, prop=0.8, group_variables=c("plot", "season", "disease"))

split_data$cal_data$id
[1] 7 19 13 25 37 14 38 20 32 26 3 21 15 9 27 10 34 16
[19] 4 40 11 35 5 17 41 12 36 6 18 42 73 55 67 49 43 44
[37] 68 56 62 50 69 75 51 45 63 52 58 82 76 70 53 71 59 77
[55] 47 54 48 66 84 72 85 109 121 91 103 86 116 98 110 92 93 123
[73] 117 99 111 112 118 94 88 100 89 113 101 95 107 102 90 120 126 114
split_data$val_data$id
[1] 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
[19] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

table(split_data$cal_data$plot,split_data$cal_data$season, split_data$cal_data$disease)
, , = 0
1 2 3 4 5 6
plot1 5 5 5 5 5 5
plot2 5 5 5 5 5 5
plot3 0 0 0 0 0 0
, , = 1
1 2 3 4 5 6
plot1 0 0 0 0 0 0
plot2 0 0 0 0 0 0
plot3 5 5 5 5 5 5

table(split_data$val_data$plot,split_data$val_data$season, split_data$val_data$disease)
, , = 1
1 2 3 4 5 6
plot3 6 6 6 6 6 6

Why not using the default method from the pls package to determine the number of components, ie the function selectNcomp? I think it would simplify the code and the manuscript. See for example p13 of the pls manual. https://cran.r-project.org/web/packages/pls/vignettes/pls-manual.pdf . But maybe you didnt like that?

To address @alistairrogers comment about examples that arent dependent on external tools, we will add a built in dataset example and associated script. Plan is to load the data using the R data() command then the rest of the steps would look the same; that is no reliance on EcoSIS

Check we use all loaded packages and all package versions are somewhere.

https://github.blog/2020-12-15-token-authentication-requirements-for-git-operations/

This does not seem to be working for me now:
devtools::source_url("https://raw.githubusercontent.com/TESTgroup-BNL/PLSR_for_plant_trait_prediction/master/R_Scripts/functions.R")

CC @JulienLamour @neo0351

In our example scripts, particularly when we create the side-by-side cal/val plots

I found a bug where the R2 doesn't match the R2 in the final validation plot

Thats because for the cal/val plotting we are putting in R2 using

val.R2 <- round(pls::R2(plsr.out,newdata=val.plsr.data)[[1]][nComps],2)

but that doesn't account for the fact that the intercept is included so instead of say giving me the results for nComp=5 its giving me the value at ncomp=4. We need to do

val.R2 <- round(pls::R2(plsr.out,newdata=val.plsr.data)[[1]][nComps+1],2)

which is what we do for the validation plot. So we need to make a change in our scripts to address this and re-run the examples :/

@julien is this also what you saw? Or was that a different issue?

@JulienLamour noticed an error in the calculation, Needs to be changed to match our reviewer responses and to match citations

I am trying to fit a PLSR model of hyperspectral (VNIR) data on SLA (mean of canopy values per tree) using psl and spectratrait and following the tutorial provided in it. But I find a lack of congruence between the results of spectratrait::find_optimal_components and what I see in the training, crossvalidation and validation results with different number of components.
find_optimal_components always select 1 to 3 components (depending on the split proportion), but all three data sets seem to improve significantly at least until 8 components.

There is surely something I am missing, but I cannot get the reason why, but there seems to be a general agreement between train and validation results, while cross-validation results have a very different pattern.
My data come from a complex experimental design with factors plot, disease and season. I use these factors for creating the validation data set...but it seems to me, that find_optimal_components does not account for this stratification. Could this be the reason for the differing behaviour of crossvalidation? Could I use spectratrait::create_data_split and then run the model iteratively for different splits to select the optimal number of components?
Even disregarding this difference, seems that the optimal number of components should be something around 8 (if we use cross-validation as a reference), but the procedure always gives values between 1 and 3 (tried with different methods).
May the lack of stratification in cross-validation be behind this weird behaviour?
Best regards,
Asier

Started but we could use more

devtools::test("PLSR_for_plant_trait_prediction")

> devtools::test("PLSR_for_plant_trait_prediction")
Loading spectratrait
Testing spectratrait
✔ |  OK F W S | Context
⠏ |   0       | Grabbing data via the EcoSIS API works                                                                                                              [1] "**** Downloading Ecosis data ****"
✔ |   1       | Grabbing data via the EcoSIS API works [8.0 s]

══ Results ════════════════════════════════
Duration: 8.0 s

[ FAIL 0 | WARN 0 | SKIP 0 | PASS 1 ]

Should we just go ahead and rename this repo to the name of the R package it builds? The current name is a holdover from when this was just a repo of scripts. It may make sense to just rename to "spectratrait" Thoughts?

@JulienLamour @regnans @neo0351 @DavidsonKen

Replace use of reshape2::melt() with new suggested package data.table::melt()

reshape2 is apparently going the way of the dodo at some point soonish, or has, or maybe no longer?? Need to confirm and either address this issue or close

Forgot to add these to the final graph

Before we finish we need to be sure both the script and the Markdown script match.

After thinking about it, we need to provide help on installing the package and depends and then only load whats needed in the scripts. Current approach is not a good one or R compliant.

I called this "How_to_PLSR" as a placeholder. We should probably change this at some point so its a bit more descriptive and matches the paper

Both approaches are less than desirable for different reasons

1. pls jackknife selection often results in more components than necessary

2. T.test permutation requires a large enough number of iterations otherwise the auto-select chooses far too few. Or if there is a lot of variance in the data it can result in the selection of a small number of components. For example with the https://ecosis.org/package/leaf-reflectance-plant-functional-gradient-ifgg-kit dataset. My tests keep selected 2 components which isnt correct. I wonder if we can refine how we find the minimum?

COuld use this to show our code works with image spectra as well

https://ecosis.org/package/canopy-spectra-to-map-foliar-functional-traits-over-neon-domains-in-eastern-united-states

General reference
https://docs.github.com/en/actions/quickstart

References for use with R code

https://usethis.r-lib.org/reference/github_actions.html
usethis::use_gh_actions_ci()

https://github.com/r-lib/actions

Make sure to create the PDF manual

in the package folder, run:

R

devtools::build_manual(pkg = ".", path = NULL)

line 117: val.pls.data<-plsr_data[!plsr_data$id %in% cal.plsr.data$id,] change to val.plsr.data<-plsr_data[!plsr_data$id %in% cal.plsr.data$id,]

line 153: object cal_hist_plot returns error. I think this is because the object cal.plsr.data which is used to make the plot is an S3 grouped object and thus qqplot doesn't know how to deal with it. I changed cal.plsr.data to a data.frame and it worked fine for the rest of the script

line 193 and 195: Error in xy.coords(x, y, setLab = FALSE) : 'x' and 'y' lengths differ. Looks like there is an additional entry "id" in the cal.plsr.data$Spectra vector which is causing the differing lengths

line 445, 449, 593, 595 and 599,: Error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ. I think this is the same issue as above. The vectors "coefs" and "vips" have a length of 1902, which is 1 longer than the length of seq(Start.wave,End.wave,1). The named last entry in both of these vectors is id

Somehow we still have

devtools::install_github(repo = "TESTgroup-BNL/PLSR_for_plant_trait_prediction", 
dependencies=TRUE)

and this should be

devtools::install_github(repo = "TESTgroup-BNL/spectratrait", dependencies=TRUE)

FYI @JulienLamour when playing with another dataset from EcoSIS I realized that your dplyr method for splitting data assumes the dataset contains "Sample_ID" which is based on the original example but not universal

    } else if (approach=="dplyr")
      cal.plsr.data <- plsr_data %>% 
        group_by_at(vars(all_of(group_variables))) %>% 
        slice(sample(1:n(), prop*n())) %>% 
        data.frame()
      val.plsr.data <- plsr_data[!plsr_data$Sample_ID %in% cal.plsr.data$Sample_ID,]

We will need to revise how we select out the not cal data so it doesnt depend on knowing what variable is in the dataset. currently I am not sure yet how do do that with that approach as it doesnt create a new variable to track which rows are part of the cal selected dataset, like the original version as revised by @neo0351 . We will need to fix this to make it more generic

Originally posted by Angie on the wrong git: serbinsh/SSerbin_etal_2019_NewPhytologist#1

Reposted here

When running the script via the RPubs version, nComps was not generated automatically (pressDFres not found in the t.test script). I had to manually set nComps <- 14 (guided by the annotation) to move to the next step.

When running the script from GitHub using RStudio, I could see the graphics provided in the Rmd file but running the code chunks didn't generate my own graphics in the graphics window of RStudio. I expected I would see the graphics there.

I tried both of these scripts from GitHub:
run_lma_plsr_example.Rmd
run_lma_plsr_example_neon.Rmd

I also think the uninitiated user would benefit from a bit more description of what the code is doing, throughout.

Warning: qplot() was deprecated in ggplot2 3.4.0.

e.g.

cal_hist_plot <- qplot(cal.plsr.data[,paste0(inVar)],geom="histogram",
                       main = paste0("Cal. Histogram for ",inVar),
                       xlab = paste0(inVar),ylab = "Count",fill=I("grey50"),col=I("black"),
                       alpha=I(.7))

Already generating the output needed to make non-parametric CIs instead of using basic pred intervals.

Just noticed some neat datasets on EcoSIS with novel measurements

e.g needle species with chemical trait data
https://ecosis.org/package/missoula-montana-lodgepole-pine---big-sagebrush-time-series

There are some wrong xlab in the coefs plots. they go from 0 to 1500 instead of 500 to 2500

@regnans
#69 (comment)

Assigning updating the README to you so you can add the suggested changes you were thinking about. Feel free to make the edits right in GH.

https://github.com/TESTgroup-BNL/PLSR_for_plant_trait_prediction/blob/Rpackagify/R/f.plot.coef.R

Cannot download data using read.csv or read_csv. Output for both provided. However if i paste the URL into Chrome, it downloads the file.

dat_raw <- read.csv(ecosis_file, na.strings = "")
Error in file(file, "rt") : 
  cannot open the connection to 'https://ecosis.org/api/package/5617da17-c925-49fb-b395-45a51291bd2d/export?metadata=true'
In addition: Warning message:
In file(file, "rt") :
  URL 'https://ecosis.org/api/package/5617da17-c925-49fb-b395-45a51291bd2d/export?metadata=true': status was 'SSL peer certificate or SSH remote key was not OK'

dat_raw <- readr::read_csv(ecosis_file)
Error in open.connection(con, "rb") : 
  SSL certificate problem: certificate has expired

@JulienLamour pointed out here that we are hard-coded a minimum number of obs needed for data split but then provide no error handling when the obs are too low. We should instead either stop and mention that the groups are insufficient or provide a solution, such as putting the obs into either cal or val randomly and then reporting as such since again this is really about selecting cal./val datasets

https://github.com/TESTgroup-BNL/spectratrait/blob/4e472c6ae0e5833faae7d2f0f2b2fd9a2b48c196/R/create_data_split.R#L30

I'm trying to be "that reviewer" that has no idea what they are doing with GitHub or R. (not that far from the truth).

Running spectra-trait_reseco_lma_plsr_example.R
I had some issues getting started, but will make that a different issue.
Once I had a the code running in R Studio (v1.4, with R 4.0.3), these are potential issues.

At Line 31: install devtools. Outputs list of packages (n = 56). Instructions about "enter one or more numbers" seem to be repeated and confusing.

…
56: readr (1.3.1 -> 1.4.0 ) [CRAN]

Enter one or more numbers, or an empty line to skip updates:list.of.packages <- c("pls","dplyr","reshape2","here","plotrix","ggplot2","gridExtra","spectratrait")
Enter one or more numbers, or an empty line to skip updates:invisible(lapply(list.of.packages, library, character.only = TRUE))
Enter one or more numbers, or an empty line to skip updates:#--------------------------------------------------------------------------------------------------#
Enter one or more numbers, or an empty line to skip updates:

Entered “1”
Next output message:

There are binary versions available but the source versions are later:
binary source needs_compilation
waldo 0.2.3 0.2.4 FALSE
promises 1.1.1 1.2.0.1 TRUE

Do you want to install from sources the package which needs compilation? (Yes/no/cancel)

Entered “Yes” (maybe guidance to this should be given in instructions).

Continued running chunks, all working until Line 99.

Error in sample_info %>% select(Plant_Species = Latin Species, Species_Code = species code, :
could not find function "%>%"

Tried to load packages that might help with "%>%"

library(dplyr)
Error: package or namespace load failed for ‘dplyr’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
namespace ‘rlang’ 0.4.7 is already loaded, but >= 0.4.9 is required

library(magrittr)
Error in value[3L] :
Package ‘magrittr’ version 1.5 cannot be unloaded:
Error in unloadNamespace(package) : namespace ‘magrittr’ is imported by ‘tibble’, ‘testthat’ so cannot be unloaded

Testing ended at this point. Unable to continue.
I have recently used dplr in other code without an issue. Not sure what has changed.

Change LMA to SLA in filename etc for expanded_spectra-trait_kit_lma_plsr_example.R
Rename files without 'expanded' and delete simple file.

After updating README

@neo0351

FYI -

> apply(plsr_data[, split_var], MARGIN = 1, FUN = function(x) paste(x, collapse = " "))
Error in apply(plsr_data[, split_var], MARGIN = 1, FUN = function(x) paste(x,  : 
  dim(X) must have a positive length

I found a bug in your split function when trying out a new dataset
https://ecosis.org/package/leaf-reflectance-plant-functional-gradient-ifgg-kit
ID: 3cf6b27e-d80e-4bc7-b214-c95506e46daa

Not yet sure what the issue is but it looks at as coded it assumes there will be two grouping variables. We need this to be flexible enough to handle 1+ grouping variables

And also FYI - if I try with two grouping vars I get this

> split_data <- create_data_split(approach=method, split_seed=2356812, prop=0.8, 
+                                 group_variables=c("Growth_Form","Plant_Species"))
NA   Cal: 79.9729364005413%
 Not enough observations

Again these functions need to be general to allow for flexibility. We will need to fix this to allow for different numbers of grouping variables with different numbers of obs.

Given the increased complexity and interest in using the example PLSR R script, the question now is whether we should create a small R package containing the functions developed for the scripts in order to make it easier to find and use these functions in custom fitting scripts? The idea would be this package would be a companion to the standard PLS package containing custom functions to make the approach easier. It could also contain 1 or 2 example datasets.