Package Functions

[lilyclements]

Introduction

This is looking at the gosset and PlackettLuce R packages so far:

From current investigation, there are different classes: Two types of ranking classes, and several types of models:

1. rankings
2. grouped_rankings
3. PlackettLuce
4. pltree
5. pladmm, bttree, gnm - to explore

Things we can do with Rankings:

Creating the Object
We create a rankings object using either rankings, rank_numeric, or rank_tricot(..., group = FALSE)

Describing the Object
For completeness sake of the functions in gosset and in the PlackettLuce packages:

adjacency(rank_obj)
aggregate(rank_obj)
choices(rank_obj)
connectivity(rank_obj)
kendallTau(rank_obj_1, rank_obj_2)
# kendallW?
set_paircomp(rank_obj)  # creates a new rank object
set_binomialfreq(set_paircomp(rank_obj))
set_binomialfreq(rank_obj)
group(rank_obj, ....) # set it to be a group ranking object!

Modelling the Object

PlackettLuce(rank_obj) # model it

Things we can do with Grouped Rankings:

Creating the Object
We create a rankings object using either group(...) to a rankings object, or rank_tricot(..., group = TRUE)

Describing the Object
For completeness sake of the functions in gosset and in the PlackettLuce packages:

kendallTau(rank_obj_1, rank_obj_2)
# kendallW?
set_paircomp(rank_obj)  # creates a new rank object
set_binomialfreq(set_paircomp(rank_obj))
set_binomialfreq(rank_obj)

Modelling the Object

PlackettLuce(rank_obj) # model it
pltree(group_rank_obj) # model it

Things we can do with PlackettLuce Models:

Creating the Object
Created by running PlackettLuce(rank_obj) to either a rankings or grouped rankings object.

Fitting the Model / Looking at the Model

worth_bar(rank_pl)
resample(rank_pl)
coef(rank_pl)
pairwise_probs(coef(rank_pl))
reliability(rank_pl, ref= "ALS 0532-6")
compare(rank_pl, rank_pl2)
anova(rank_pl)
AIC(rank_pl)
deviance(rank_pl)
summary(rank_pl)
fitted(rank_pl)

# then unsure about:
# kendallTau? kendallW?
#vcov(rank_pl)
#qvcalc(rank_pl)
#simulate(rank_pl)
#PlackettLuce::itempar(rank_pl, rank_pl)
#crossvalidation(rank_pl ~ ., data = R)
#akaike_weights(rank_pl)
#pseudo_rank(rank_pl)

Comparing Models

# Let rank_pl_list be a list of PL models:
worth_map(rank_pl_list)

# Then let models be a vector of AICs:
models = c(mod1 = AIC(rank_pl), mod2 = AIC(rank_pl2))
akaike_weights(models)

Things we can do with pltrees Models:

Creating the Object
Created by running pltree(rank_obj) to a grouped rankings object.

Fitting the Model / Looking at the Model

reliability(mod, ref = "A")
plot(mod)
#crossvalidation(mod, data = R)
AIC(mod)
deviance(mod)
worth_map(mod)
worth_bar(mod)
regret(mod)
coef(mod)
predict(mod)
fitted(mod)
itempar(mod)
#qvcalc(mod)
summary(mod)
kendallTau(R, predict(mod))
kendallW(R, predict(mod))
pairwise_probs(coef(mod))

Comparing Models

# Let list_mods be a list of pltree models:
list_mods <- list(mod, mod2)
worth_map(list_mods, labels = c("1", "2"))

# Then let models be a vector of AICs:
models = c(mod1 = AIC(rank_pl), mod2 = AIC(rank_pl2))
akaike_weights(models)

[lilyclements]

For the chirps package, is there anything in there that we want that we don't have yet?

We want to look at the ClimMob-Analysis repo

[lilyclements]

In practice, CHIRPS would be the preferred method in their training. We want to investigate CHIRPS - what do we have, what do we want that we do not have currently in R-Instat?

[lilyclements]

ag5Tools has two functions:

1. ag5_download
2. ag5_extract

Are these used in Tricot training? Should these be included in the R-Instat menu?

[lilyclements]

They are not used in training - they require an account, involve data treatment, and could potentially open up a lot of complications that may distract from the core workflow. But, they could be useful in real-wrld scenarios so might be sensible to include

[lilyclements]

The ClimMobTools Package:

1. There are a set of functions useful in importing from ClimMob, with a related issue here.
2. The getProjectProgress function is deprecated and to be removed.
3. Then there is the randomize function, which might be relevant to the discussion here so is added. This randomisation method is the core of the whole tricot approach. This is used to illustrate how randomisation is done – just show the result is of interest. Not sure how to incorporate it, unless it is as one of Roger's examples in the help file area.
4. rankTricot is the same as rank_tricot (see R code below, and confirmed).
5. rmGeoIdentity - gives error: could not find function "st_point". This is a personal function for the pkg author, not used in the training. It removes the identity, replaced co-ord with boundary not an exact value. I'm not sure if that is something we already have in R-Instat or would be interested in, so I'll just flag it up here.
6. getTraitList is used in the vignettes, but there is a question on if that is needed in R-Instat if we do a "Define" dialog

library(ClimMobTools)
library(gosset)

data("beans", package = "PlackettLuce")


R <- rankTricot(data = beans,
                items = c(1:3),
                input = c(4:5))
R1 <- rank_tricot(beans,
                 items = c(1:3),
                 input = c(4:5), 
                 group = FALSE)
identical(R, R1)


G = rankTricot(data = beans,
               items = c(1:3),
               input = c(4:5),
               group = TRUE,
               additional.rank = beans[c(6:8)])
G1 = rank_tricot(data = beans,
                 items = c(1:3),
                 input = c(4:5),
                 group = TRUE,
                 additional.rank = beans[c(6:8)])
identical(G, G1)

[lilyclements]

TODO: look into nasapower package

[lilyclements]

In training, they use nasapower as it has a simpler workflow than ag5Tools and CHIRPS. I suggest then we do want to look into this nasapower package and have it linked into for importing.
We should also have ag5Tools in there too since they would use that outside of training
In practice, CHIRPS would be preferred to nasapower. We should investigate CHIRPS: what’s available to us, and what we may want to integrate going forward.

[lilyclements]

The purpose of the climatrends package is to bridge from the climate data to Climmob (the tricot data).

In training, they typically fit four models to demonstrate how results can change depending on different factors. The idea is to emphasise that data has context (OxC), and results can vary based on model choices:

a) No covariates (to show heterogeneity and how rankings change with factors).
b) Demographic variables (which lead to another set of results).
c) A climatology, such as rainfall, for a third result.
d) A model using demographic and climatology covariates to get a better result

In the third model, climatrends come into play. The functions they use for this are GDD, rainfall, and temperature.
I note that @rdstern mentioned that "the rainfall indices are taken from the climdex suite, so they are already in R-Instat. The temperature indices are really interesting. As we proceed I would guess we should be adding them into R-Instat, rather than adding the package.". So adding the temperature indices in that way seems sensible, as was Roger's thoughts before the meeting!

Other functions in climatrends (like cumdrought, cumrain, and crop_sensitive) are not a focus in tricot, though important for other cases. To what extent do we want to dig into these?

[rdstern]

@jkmusyoka should also be looking at the chirps package for our climatic work in general and particularly for our uptfilling work. So far, in our outfilling he has supplied the chirps (and tamsat, etc) data in a file ready to add to the station data. Maybe the chirps package can make the use of chirps easier.

However, what I think this all implies for the tricot work is that they currently largely use chirps, simply as it stands. I suggest Danny's paper shows that chirps (and other satellite estimates) may not be great for the sort of indices, needed by PICSA - and I suggest they will be similar to useful indices needed by tricots?

Now Emily's work proposes a locally derived adjustment to chirps, so it is more appropriate for that province or country. It is then used to outfill, but can equally be used to adjust chirps when there are no ground stations. And ground stations don't seem to be mentioned by the TRICOTS people.

So, maybe they could be excited by the results of this work, that could improve the climatic indices they use. If so - and after discussions with Emily, we may keep this in reserve for a bit, but it could become a major enhancement of "our" work, and probably not part of this initial contract.

[lilyclements]

Very interesting! Thanks for looping me in on how this is all fitting together.

I'll leave the chirps package up to @jkmusyoka - do let me know your findings on this. Sounds very interesting!

[jkmusyoka]

@lilyclements @rdstern
R-Instat already has a dialog for importing chirps data. It seems to be working alright though it takes time to run for large downloads and often times out.

Looking at the code, R-Instat seems to read the chirps data from the IRI website which I believe also refers to the chirps website for the estimates.
data_book$download_from_IRI(source="UCSB_CHIRPS", data="daily_improved_global_0p05_prcp", path="C:/Users/jkmus/OneDrive - IDEMS International/Documents", min_lon=34.58, max_lon=35.08, min_lat=0, max_lat=0.5, min_date=as.Date(x="2021/1/1"), max_date=as.Date(x="2021/1/31"), name="ucsb chirps", download_type="Point")

Unlike chirps package, R-Instat's implementation returns a .T column for date (I believe its special date-time column thing perhaps Julian?). The Make Data dialog is then needed to create a proper date column from this .T column. So this could be an improvement we may want to make with the current implementation.

[lilyclements]

@jkmusyoka great, thank you for this. They only use chirps for importing from what I can see, and perhaps chirts as well.

1. Do we have a way to get chirts data in?

2. I don't think we have the flexibility to do what they're doing here - i.e., import for different longitudes and latitudes, and for different planting dates, that are read from a dataset:

library(chirps)
library(gosset)

data(nicabean)

# I'm just going to show this example for the first six rows of covar, because this takes about an hour to run if we do it for all the rows!
# But in covar we have planting_date, latitude, and longitude. 
# This is used for getting the chirps data for different longitudes/latitudes/planting dates:

covar <- head(nicabean$covar)

dates = c(min(covar[, "planting_date"]),
          max(covar[, "planting_date"]) + 70)

chirps = get_chirps(covar[, c("longitude","latitude")], 
                    dates = as.character(dates),
                    as.matrix = TRUE,
                    server = "ClimateSERV")
# rename the date variable
newnames = dimnames(chirps)[[2]]
newnames = gsub("chirps-v2.0.", "", newnames)
newnames = gsub("[.]", "-", newnames)
dimnames(chirps)[[2]] = newnames

Can you see a way we can do this in R-Instat?

[jkmusyoka]

@lilyclements my responses to your two points as follows:-

1. R-Instat does not yet have a dialog for importing chirts. Only rainfall. If needed we may need build this dialog from scratch
2. I have found a way of using R-Instat dialogs to download chirps data values for multiple stations, long and lat just like in the R script.
   Step 1 - The starting point is the 2015 chirps dataset. This netCDF file must be available in your computer. The file can be downloaded for a region of interest using an R script. R-Instat does not have a way of downloading this file at the moment
   Step 2 - In R-Instat, go to Climatic>File>Import and Tidy NetCDF. Browse to the netCDF file in your computer
   
   Step 3 - Click on Options in the dialog to open a sub-dialog box. Select Points and supply the lon, lat and date columns from the covar file already open in R-Instat.
   
   Step 4 - Click Return and then OK to complete the extraction.
   A dataframe is returned with the Chirps values in the precip column
   

I am yet to explore the process for dates spanning more than 1 year. Should I be checking and if yes do you have an example R script?