data("userlogs")data.userlogs.raw= userlogshead(data.userlogs.raw)

Average task completion time (arithmetic mean) for eachcondition.

We need a few helper functions to so that they take the samearguments and return the same output.These functions will help uscalculate the mean point estimate and the upper and lower bounds of the95% confidence interval using the bootstrap method and t-testmethod.

bootstrappedCI<-function(observations, conf.level,seed =0) {  samplemean<-function(x, d) {return(mean(x[d]))}  pointEstimate<-samplemean(observations)if (!(is.na(seed)|is.null(seed))){set.seed(seed)# make deterministic  }  bootstrap_samples<- boot::boot(data = observations,statistic = samplemean,R =5000)  bootci<- boot::boot.ci(bootstrap_samples,type ="bca",conf = conf.level)c(pointEstimate,  bootci$bca[4], bootci$bca[5])}tCI<-function(observations, conf.level) {  pointEstimate<-mean(observations)  sampleSD<-sd(observations)  sampleN<-length(observations)  sampleError<-qt(1-(1-conf.level)/2,df = sampleN-1)* sampleSD/sqrt(sampleN)c(pointEstimate, pointEstimate- sampleError, pointEstimate+ sampleError)}

Next we initialise the multiverse object within which this analysiswill take place.

M=multiverse()

Now we define the parameters we want to consider in the multiverse:confidence interval method (ci_method) and datatransformation method (data_transform). We also define aparameter for confidence level, as the choice of a 95% confidence levelis arbitrary and we can choose to instead present our results withalternate confidence levels. Here we vary between: 50%, 89%, 95%, and99.9%. Thus our multiverse consists of: 2\(\times\) 2\(\times\) 7 different analysiscombinations.

```{multiverse default-m-1, inside = M}ci_method <- branch(ci_method,    "t based"   ~ "tCI",    "bootstrap" ~ "bootstrappedCI")data_transform <- branch(data_transform,    "log-transformed" ~ log,    "untransformed" ~ identity)conf_level <-  branch(conf_level,    "50%" ~ 0.5,    "89" ~ 0.89,    "95%" ~ 0.95,    "99%" ~ 0.99)```

expand(M)

## # A tibble: 16 × 8##    .universe ci_method data_transform  conf_level .parameter_assignment##        <int> <chr>     <chr>           <chr>      <list>               ##  1         1 t based   log-transformed 50%        <named list [3]>     ##  2         2 t based   log-transformed 67         <named list [3]>     ##  3         3 t based   log-transformed 95%        <named list [3]>     ##  4         4 t based   log-transformed 99%        <named list [3]>     ##  5         5 t based   untransformed   50%        <named list [3]>     ##  6         6 t based   untransformed   67         <named list [3]>     ##  7         7 t based   untransformed   95%        <named list [3]>     ##  8         8 t based   untransformed   99%        <named list [3]>     ##  9         9 bootstrap log-transformed 50%        <named list [3]>     ## 10        10 bootstrap log-transformed 67         <named list [3]>     ## 11        11 bootstrap log-transformed 95%        <named list [3]>     ## 12        12 bootstrap log-transformed 99%        <named list [3]>     ## 13        13 bootstrap untransformed   50%        <named list [3]>     ## 14        14 bootstrap untransformed   67         <named list [3]>     ## 15        15 bootstrap untransformed   95%        <named list [3]>     ## 16        16 bootstrap untransformed   99%        <named list [3]>     ## # ℹ 3 more variables: .code <list>, .results <list>, .errors <list>

```{multiverse default-m-2, inside = M}duration <- do.call(data_transform, list(data.userlogs.raw$duration))```

```{multiverse default-m-3, inside = M}modality <- data.userlogs.raw$modalitynameci.physical_notouch <- do.call(ci_method, list(duration[modality == 'physical-notouch'], conf_level))ci.physical_notouch <- setNames(as.list(c("physical_notouch", ci.physical_notouch)), c("modality", "estimate", "conf.low", "conf.high"))ci.physical_touch <- do.call(ci_method, list(duration[modality == 'physical-touch'], conf_level))ci.physical_touch <- setNames(as.list(c("physical_touch", ci.physical_touch)), c("modality", "estimate", "conf.low", "conf.high"))ci.virtual_prop <- do.call(ci_method, list(duration[modality == 'virtual-prop'], conf_level))ci.virtual_prop <- setNames(as.list(c("virtual_prop", ci.virtual_prop)), c("modality", "estimate", "conf.low", "conf.high"))ci.virtual_mouse <- do.call(ci_method, list(duration[modality == 'virtual-mouse'], conf_level))ci.virtual_mouse <- setNames(as.list(c("virtual_mouse", ci.virtual_mouse)), c("modality", "estimate", "conf.low", "conf.high"))df <- rbind.data.frame(ci.physical_notouch, ci.physical_touch, ci.virtual_prop, ci.virtual_mouse, make.row.names = FALSE, stringsAsFactors = FALSE)df <- transform(df, estimate = as.numeric(estimate), conf.low = as.numeric(conf.low), conf.high = as.numeric(conf.high))```

execute_multiverse(M)

Extracting the results from the multiverse object

Next, we need to extract the results from the multiverse. The resultsfor each unique analysis combination (a universe in our multiverse), isstored in an environment in the.results column. We canextract data frames from this column using the function. This creates anew column in our data frame,summary which itself consistsof data frames.

df.mtbl<-expand(M)df.mtbl$summary=map(df.mtbl$.results,"df")head(df.mtbl)

## # A tibble: 6 × 9##   .universe ci_method data_transform  conf_level .parameter_assignment##       <int> <chr>     <chr>           <chr>      <list>               ## 1         1 t based   log-transformed 50%        <named list [3]>     ## 2         2 t based   log-transformed 67         <named list [3]>     ## 3         3 t based   log-transformed 95%        <named list [3]>     ## 4         4 t based   log-transformed 99%        <named list [3]>     ## 5         5 t based   untransformed   50%        <named list [3]>     ## 6         6 t based   untransformed   67         <named list [3]>     ## # ℹ 4 more variables: .code <list>, .results <list>, .errors <list>,## #   summary <list>

As we can see above, each row in thesummary columnconsists of a\(4 \times 4\), dataframe which we will need to unpack. We will use the function to expandthe different columns of the data frame into their own columns. Finally,we use the function to unnest the rows of the data frame into their ownrow. Note that we have a.universe column which indexeseach universe in our multiverse i.e. each unique analysiscombination.

Below we can see the result of this transformation. You can see thatwe have created four new columns (modality,estimate,conf.low,conf.high).In addition, we have four rows for each universe corresponding to theresults for each of the four conditions in our experiment.

df.mtbl<-unnest_wider(df.mtbl,c(summary))df.mtbl<-unnest(df.mtbl,cols =c(modality, estimate, conf.low, conf.high))head(df.mtbl)

## # A tibble: 6 × 12##   .universe ci_method data_transform  conf_level .parameter_assignment##       <int> <chr>     <chr>           <chr>      <list>               ## 1         1 t based   log-transformed 50%        <named list [3]>     ## 2         1 t based   log-transformed 50%        <named list [3]>     ## 3         1 t based   log-transformed 50%        <named list [3]>     ## 4         1 t based   log-transformed 50%        <named list [3]>     ## 5         2 t based   log-transformed 67         <named list [3]>     ## 6         2 t based   log-transformed 67         <named list [3]>     ## # ℹ 7 more variables: .code <list>, .results <list>, .errors <list>,## #   modality <chr>, estimate <dbl>, conf.low <dbl>, conf.high <dbl>

Visualising the results

We will then sort the results, and transform the log transformedvariables back on to the natural scale. We are then ready to visualisethe result by animating over each universe.

df.mtbl<-arrange(df.mtbl, conf_level,desc(data_transform),desc(ci_method))df.results<- df.mtbldf.results$estimate[df.mtbl$data_transform=="log-transformed"]=exp(df.mtbl$estimate[df.mtbl$data_transform=="log-transformed"])df.results$conf.low[df.mtbl$data_transform=="log-transformed"]=exp(df.mtbl$conf.low[df.mtbl$data_transform=="log-transformed"])df.results$conf.high[df.mtbl$data_transform=="log-transformed"]=exp(df.mtbl$conf.high[df.mtbl$data_transform=="log-transformed"])df.results|>head()

## # A tibble: 6 × 12##   .universe ci_method data_transform conf_level .parameter_assignment##       <int> <chr>     <chr>          <chr>      <list>               ## 1         5 t based   untransformed  50%        <named list [3]>     ## 2         5 t based   untransformed  50%        <named list [3]>     ## 3         5 t based   untransformed  50%        <named list [3]>     ## 4         5 t based   untransformed  50%        <named list [3]>     ## 5        13 bootstrap untransformed  50%        <named list [3]>     ## 6        13 bootstrap untransformed  50%        <named list [3]>     ## # ℹ 7 more variables: .code <list>, .results <list>, .errors <list>,## #   modality <chr>, estimate <dbl>, conf.low <dbl>, conf.high <dbl>

p<- df.results|>ggplot()+geom_vline(xintercept =0,colour ='#979797' )+geom_point(aes(x = estimate,y = modality))+geom_errorbarh(aes(xmin = conf.low,xmax = conf.high,y = modality),height =0)+transition_manual( .universe )+theme_minimal()animate(p,nframes =28,fps =2)

## `nframes` and `fps` adjusted to match transition

The figure above shows the (geometric) mean completion time for eachcondition. At first sight, physical touch appears to be consistentlyfaster than the other conditions, across all possible analysiscombinations specified in the multiverse. However, since condition is awithin-subject factor, it is preferable to examine within-subjectdifferences, which we show in the next section.

Sidebar: Using the tidyverse to extract and visualize theresults

So far we have tried to keep this analysis in base R. However, thesteps involved in extracting and visualising results from the multiversemaybe more convenient for some using the tidyverse API. We can implementthe steps that we have taken in the previous three code block in thefollowing way using the tidyverse. The result of this is the same dataframe that was created in the previous step, and was used as input datato ggplot,df.results. Themultiverse packageis flexible, and does not restrict you towards a particular style of Rprogramming.

expand(M)|>mutate(summary =map(.results,"df") )|>unnest_wider(c(summary))|>unnest(cols =c(modality, estimate, conf.low, conf.high))|>mutate(estimate =ifelse(data_transform=="log-transformed",exp(estimate), estimate),conf.low =ifelse(data_transform=="log-transformed",exp(conf.low), conf.low),conf.high =ifelse(data_transform=="log-transformed",exp(conf.high), conf.high)  )|>arrange(conf_level,desc(data_transform),desc(ci_method))

## # A tibble: 64 × 12##    .universe ci_method data_transform  conf_level .parameter_assignment##        <int> <chr>     <chr>           <chr>      <list>               ##  1         5 t based   untransformed   50%        <named list [3]>     ##  2         5 t based   untransformed   50%        <named list [3]>     ##  3         5 t based   untransformed   50%        <named list [3]>     ##  4         5 t based   untransformed   50%        <named list [3]>     ##  5        13 bootstrap untransformed   50%        <named list [3]>     ##  6        13 bootstrap untransformed   50%        <named list [3]>     ##  7        13 bootstrap untransformed   50%        <named list [3]>     ##  8        13 bootstrap untransformed   50%        <named list [3]>     ##  9         1 t based   log-transformed 50%        <named list [3]>     ## 10         1 t based   log-transformed 50%        <named list [3]>     ## # ℹ 54 more rows## # ℹ 7 more variables: .code <list>, .results <list>, .errors <list>,## #   modality <chr>, estimate <dbl>, conf.low <dbl>, conf.high <dbl>

execute_multiverse(M)

Extracting and visualizing the results from the multiverse

We then use the workflow described in the previous section to extractresults for each universe in the multiverse. We then use gganimate toplot the results.

df.results.diff<-expand(M)|>extract_variables(df.diffs)|>unnest(c(df.diffs))|>arrange(desc(data_transform), conf_level,desc(ci_method))df.results.diff|>head()

## # A tibble: 6 × 12##   .universe ci_method data_transform conf_level .parameter_assignment##       <int> <chr>     <chr>          <chr>      <list>               ## 1         5 t based   untransformed  50%        <named list [3]>     ## 2         5 t based   untransformed  50%        <named list [3]>     ## 3         5 t based   untransformed  50%        <named list [3]>     ## 4        13 bootstrap untransformed  50%        <named list [3]>     ## 5        13 bootstrap untransformed  50%        <named list [3]>     ## 6        13 bootstrap untransformed  50%        <named list [3]>     ## # ℹ 7 more variables: .code <list>, .results <list>, .errors <list>,## #   modality <chr>, estimate <dbl>, conf.low <dbl>, conf.high <dbl>

p<- df.results.diff|>ggplot()+geom_vline(xintercept =0,colour ='#979797' )+geom_point(aes(x = estimate,y = modality))+geom_errorbarh(aes(xmin = conf.low,xmax = conf.high,y = modality),height =0)+transition_manual( .universe )+theme_minimal()animate(p,nframes =28,fps =4)