In stephens999/dscr: Dynamic Statistical Comparisons in R

knitr::opts_chunk$set(collapse = TRUE,comment = "#")

First we load the dscr package, as well as the psych package used in the example below.

library(dscr)
library(psych)

This vignette is designed to illustrate the (in development) dscr package for Dynamic Statistical Comparisons in R.

We'll do a simulation study to assess methods for estimating the mean of a distribution. The simulations will be performed under three scenarios: a normal distribution, a uniform distribution and a Cauchy (t distribution on 1df). And we'll compare three methods: sample mean, sample median and the Winsorized mean (from the psych package).

First load the library, and initialize a new "dsc" using new_dsc. This object will be used to store the details of the dsc.

library(dscr)
dsc_osl = new_dsc("one_sample_location","one_samp_loc")

Now define the function to create data, which consists of input, and meta-data. The input is what the methods will be given. The meta data will be used when scoring the methods. So here, the input is a random sample, and the meta-data is the true mean (0).

This function can take a single parameter args, which is a list. In this case we use args to pass number of samples and the type distribution to be simulated from.

datamaker = function(args){

  nsamp=args$nsamp
  disttype=args$disttype

  if(disttype=="normal"){
    input = list(x=rnorm(nsamp,0,1))
    meta =  list(truemean=0)
  }

  if(disttype=="uniform"){
    input = list(x=runif(nsamp,-1,1))
    meta = list(truemean=0)
  }

  if(disttype=="Cauchy"){
    input = list(x=rt(nsamp,df=1))
    meta = list(truemean=0)
  }

  data = list(meta=meta,input=input)

  return(data)
}

Now define the scenarios that use this datamaker. Each scenario is determined by the datamaker, its arguments, and the seeds it uses. (Maybe the list of seeds should not be part of the scenario definition - I'm not sure). Each scenario also has a name which is used for naming the scenario in results and also output directories - so use something that is a valid filename.

add_scenario(dsc_osl,name="normal",datamaker,
             args=list(disttype="normal",nsamp=1000),seed=1:100)
add_scenario(dsc_osl,name="uniform",datamaker,
             args=list(disttype="uniform",nsamp=1000),seed=1:100)
add_scenario(dsc_osl,name="Cauchy",datamaker,
             args=list(disttype="Cauchy",nsamp=1000),seed=1:100)

Now define the methods. They have to have the form where they take "input" and produce "output" in a specified format. In this case the input format is a list with one component (x). The output format is a list with one component (meanest), the estimated mean.

Effectively we have to write a "wrapper" function for each of our three methods that makes sure that they conform to this input-output requirement. (Note that the winsor.wrapper function makes use of the function winsor.mean from the psych package.) Note that we allow for additional arguments ot each function, but don't use them here.

mean.wrapper = function(input,args){
  return(list(meanest = mean(input$x)))
}

median.wrapper = function(input,args){
  return(list(meanest = median(input$x)))
}

winsor.wrapper = function(input,args){
  return(list(meanest = psych::winsor.mean(input$x,trim=0.2)))
}

Now define a list of the methods we'll use. Each method is defined by its name, the function used to implement it, and any additional arguments (none here):

add_method(dsc_osl,name="mean",fn =mean.wrapper,args=NULL)
add_method(dsc_osl,name="median",fn=median.wrapper,args=NULL)
add_method(dsc_osl,name="winsor",fn=winsor.wrapper,args=NULL)

And define a score function that says how well a method has done. Here we'll use squared error and absolute error:

score = function(data, output){
  return(list(squared_error = (data$meta$truemean-output$meanest)^2,
         abs_error = abs(data$meta$truemean-output$meanest)))
}
add_score(dsc_osl,score)

Now we'll run all the methods on all the scenarios:

res=run_dsc(dsc_osl)

This returns a dataframe with the results of running all the methods on all the scenarios:

head(res)

And we can summarize the results (eg mean squared error) using the aggregate function:

aggregate(abs_error~method+scenario,res,mean)
aggregate(squared_error~method+scenario,res,mean)

Now suppose we are coming in and want to add a method, say the trimmed mean, to the comparison. Suppose also we want to try out the trimmed mean with two different settings of the trim argument. Here is what we do (note that the different settings of the argument are treated as different methods, but the two methods use the same function).

trimmedmean.wrapper = function(input,args){
  return(list(meanest=mean(input$x,trim=args$trim)))
}

add_method(dsc_osl,name="trimmedmean1",fn = trimmedmean.wrapper,
           args=list(trim=0.2))
add_method(dsc_osl,name="trimmedmean2",fn = trimmedmean.wrapper,
           args=list(trim=0.4))
res=run_dsc(dsc_osl)

aggregate(abs_error ~ method + scenario,res,mean)
aggregate(squared_error ~ method + scenario,res,mean)

Note that at present run_dsc does not recreate any files that are already there. Thus in this case it is only running the new methods (trimmedmean1 and trimmedmean2)- the results for other methods are already there. If you want to force it to recreate the files then you need to delete them manually before running run_dsc. You can do this from with R using the functions reset_dsc() to remove all the results, or reset_scenario or reset_method to remove results only for a particular method or scenario.

I will note that it is possible within run_dsc to specify to run the dsc on only a subset of the methods and scenarios.

res = run_dsc(dsc_osl,c("Cauchy","normal"),c("trimmedmean1"))

aggregate(abs_error ~ method + scenario,res,mean)
aggregate(squared_error ~ method + scenario,res,mean)

Finally, dscr has the ability to add a method that can "cheat" by taking not only the input data, but also the meta data. Sometimes this can be useful to allow implementation of a "gold standard" method against which other methods can be compared (even if the gold standard is not achievable in practice). In this case the method "truth" trivially cheats by passing back the truth as its estimate! This is of course just for illustration - in practice the gold standard method would usually be more subtle than this. For instance a gold standard method might perform (in some sense optimal) inference under the model used to simulate the data...

truth.wrapper = function(input,meta,args){
  return(list(meanest = meta$truemean))
}
add_method(dsc_osl,"truth",truth.wrapper,gold_flag=TRUE)

Then we can re-run the dsc:

res = run_dsc(dsc_osl)

aggregate(abs_error ~ method + scenario,res,mean)
aggregate(squared_error ~ method + scenario,res,mean)

stephens999/dscr documentation built on May 30, 2019, 3:20 p.m.