Appendix

censoringCode.Rmd

10 Dec 2020

## Initialization

rm(list = ls())

library("knitr")

library("survival")

library("OptimalCutpoints")

library("xtable")

## Survival functions

This function takes 2 quantile times and the corresponding quantile
proportions and returns the parameters of a Weibull distribution

@param Q_A A positive number

@param Q_B A positive number

@param p_A A probability $p_{A}\in(0,1)$ corresponding to the quantile
indicated by the time Q_A

@param p_B A probability $p_{B}\in(0,1)$ corresponding to the quantile
indicated by the time Q_B

@return a list The parameters _shape_ and _scale_ of a Weibull distribution

quantiles2weibull <- function(Q_A, Q_B, p_A = 0.5, p_B = 0.99) {

k = log(log(1 - p_B) / log(1- p_A)) / log(Q_B / Q_A)

b = Q_A / (-log(1- p_A))^(1/k)

list(shape = k, scale = b)

}

This function takes 2 quantile times and the corresponding quantile
proportions and returns a function that transform a vector of probabilities
into a vector of times following the corresponding Weibull distribution

@param Q_A A positive number

@param Q_B A positive number

@param p_A A probability $p_{A}\in(0,1)$ corresponding to the quantile
indicated by the time Q_A

@param p_B A probability $p_{B}\in(0,1)$ corresponding to the quantile
indicated by the time Q_B

@param p A vector of probabilities

@return a vector The vector of times

inverseWeibull <- function(Q_A, Q_B, p_A = 0.5, p_B = 0.99) {

function(p = runif(1)) {

qWeibull <- quantiles2weibull(Q_A, Q_B, p_A, p_B)

with(qWeibull, sapply(p, function(p) scale * (- log(p))^(1/shape)))

}

}

This function accepts a number of cases and the parameters of a mixed model
(median, $Q_{99}$, and cure-rate), estimates the cured cases and the inverse
Weibull event times, and returns a virtual dataset with complete follow-up to
be used by the survival package.

@param nCases The number of cases in the virtual dataset

@param Q_A A positive number, the median event time

@param Q_B A positive number, the $Q_{99}$ event time

@param pCure The cure rate $pCure\in(0,1)$ of the mixed model

@param x The label for the dataset

@return the virtual dataset. A data.frame with columns c(time, status, x)

completeFA <- function(nCases = 10, Q_A = 25, Q_B = 100, pCure = 0, x =
"Uncensored") {

nEvents <- as.integer(nCases * (1 - pCure))

time = sort(inverseWeibull(Q_A * (1 - pCure), Q_B)(runif(nEvents)))

tmax = max(time)

time = c(time, rep(tmax, nCases - nEvents))

status <- c(rep(1, nEvents), rep(0, nCases - nEvents))

data.frame(time, status, x = x, stringsAsFactors = FALSE)

}

## Censoring functions

This function takes a complete follow-up dataset and the parameters of a mixed
model (events median and $Q_{99}$, and proportion of censoring), estimates a
vector of random probabilities and the inverse Weibull time, updates the
dataset time and status when the new time is shorter, and returns a censored
dataset

@param uncensored The complete follow-up dataset

@param Q_A The median time of the distribution

@param Q_B the $Q_{99}$ time of the distribution

@param pCensoring A proportion of censored cases

@return a censored dataset. A data.frame with columns c(time, status, x)

timeCensoring <- function(uncensored, Q_A = 25, Q_B = 100, pCensoring = 0) {

censored <- uncensored

n <- nrow(censored)

pRandom <- runif(n)^((1-pCensoring)/(pCensoring))

timeCens <- sapply(pRandom, inverseWeibull(Q_A, Q_B))

censor <- timeCens < censored$time

censored$time[censor] <- timeCens[censor]

censored$status[censor] <- 0

censored$x <- paste("Time censoring", round(100 * pCensoring, 1), "%")

censored[order(censored$time), ]

}

This function takes a complete follow-up dataset and the parameters of a mixed
model (events median and $Q_{99}$, and proportion of censoring), estimates a
vector of random interim times (from a vector of random probabilities and the
inverse Weibull time) and of random recruit times, updates the dataset time
and status when the difference of recruit and interim times is shorter than
the original time, and returns a censored dataset

@param uncensored The complete follow-up dataset

@param Q_A The median time of the distribution

@param Q_B the $Q_{99}$ time of the distribution

@param pCensoring A proportion of censored cases

@return a censored dataset. A data.frame with columns c(time, status, x)

interimCensoring <- function(uncensored, Q_A = 25, Q_B = 100, pCensoring = 0)
{

censored <- uncensored

n <- nrow(censored)

interimTime <- inverseWeibull(Q_A, Q_B)(pCensoring) * 2

recruitTime <- runif(n, 0, interimTime)

interval <- interimTime - recruitTime

censor <- (interimTime - recruitTime) < censored$time

censored$status[censor] <- 0

censored$time[censor] <- interval[censor]

censored$x <- paste("Interim at t =", round(interimTime, 1))

censored <- censored[censored$time > 0, ]

censored[order(censored$time), ]

}

This function takes a complete follow-up dataset and a proportion of
censoring, selects a random sample of cases with that probability, updating
the status and times (shortening them a random amount), and returns a censored
dataset

@param uncensored The complete follow-up dataset

@param pCensoring A proportion of censored cases

@return a censored dataset. A data.frame with columns c(time, status, x)

caseCensoring <- function(group0, pCensoring = 0) {

n <- nrow(group0)

nCens <- as.integer(n * pCensoring)

censor <- sample(1:n, nCens)

censored <- group0

nCensored <- as.integer(nrow(group0) * pCensoring)

censored$time[censor] <- group0$time[censor] * runif(nCensored, 0.2, 1)

censored$status[censor] <- 0

censored$x <- paste("Case", round(100 * (1 - pCensoring), 1), "% censored")

censored[order(censored$time), ]

}

## Bias indexes (BI)

This function estimates the quantile bias index (QBI) of a dataset

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return a number The QBI

QBI <- function(dataset, time = "time", status = "status", event = 1, ...) {

events <- dataset[dataset[, status] == event, time]

censor <- dataset[dataset[, status] != event, time]

aCensor <- censor[censor < max(events)]

if(length(aCensor) == 0) NA else quantile(events, 0.95) / quantile(aCensor,
0.95)

}

This function estimates the scaled quantile bias index (SQBI) of a dataset

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return a number The SQBI

SQBI <- function(dataset, time = "time", status = "status", event = 1, ...) {

QBI(dataset, time, status, event) / 1.2

}

This function estimates the under the mean bias index (UMBI) of a dataset

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return a number The UMBI

UMBI <- function(dataset, time = "time", status = "status", event = 1, ...) {

events <- dataset[dataset[, status] == event, time]

censor <- dataset[dataset[, status] != event, time]

lastEvent <- max(events)

meanEvent <- mean(events)

aCensor <- censor[censor < lastEvent]

sum(aCensor < meanEvent) / length(aCensor)

}

This function estimates the adjusted bias index (ABI) of a dataset

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return a number The ABI

ABI <- function(dataset, time = "time", status = "status", event = 1, ...) {

events <- dataset[dataset[, status] == event, time]

censor <- dataset[dataset[, status] != event, time]

lastEvent <- max(events)

meanEvent <- mean(events)

aCensor<- censor[censor < lastEvent]

formula1 <- Surv(time, status == event) ~ 1

survfit1 <- survfit(formula1, data = dataset)

pLongTerm <- summary(survfit1, times = lastEvent)$surv

medianCorrection <- mean(events) / median(events)

longTermCorrection <- exp(pLongTerm)

sum(aCensor < meanEvent) / length(aCensor) * medianCorrection *
longTermCorrection

}

This function estimates the scaled adjusted bias index (ABI) of a dataset

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return

SABI <- function(dataset, time = "time", status = "status", event = 1, ...) {

ABI(dataset, time, status, event) / 0.932

}

## Experiment functions

This function extracts the hazard ratio and p-value between the two groups
labeled in the x column of a survival dataset

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@param x Labels for the two groups in the dataset

@return a list with the hazard ratio and p-value betwen the two groups

getHR <- function(dataset, time = "time", status = "status", event = 1, x =
"x") {

formula1 = Surv(time, status == event) ~ x

cx <- summary(suppressWarnings(coxph(formula1, data = dataset)))

list(hr = cx$conf.int[[1]], p = cx$coefficients[[5]])

}

This function takes a complete follow-up and a censored datasets, the names of
columns corresponding to time and status, and the label of status
corresponding to censoring, and a BI function, returning a vector with the
parameters of the trial

@param group0

@param group1

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return a vector The parameters of the trial: number of cases, proportion of
long-term survivors, hazard ratio, p-value and BI result

trialResults <-

function(group0, group1, time = "time", status = "status", event = 1) {

events <- group1[group1[, status] == event, time]

lastEvent <- max(events)

aCensor <- group1[group1[, time] < lastEvent, ]

pCensor <- sum(aCensor[, status] != event) / nrow(aCensor)

pLongTerm <- 1 - nrow(aCensor) / nrow(group1)

trialDataset <- rbind(group1, group0)

trialDataset$x <- factor(trialDataset$x, levels = c(group0$x[[1]],
group1$x[[1]]))

cox <- getHR(trialDataset)

c(nrow(group1), pCensor, pLongTerm, cox$hr, cox$p,

SQBI(group1), UMBI(group1), SABI(group1))

}

This function simulates a number of virtual trials with the indicated
conditions and returns a data.frame with the result of the trial. Every row is
a virtual trial and columns are type of censoring, number of cases, proportion
of censoring, proportion of long-term survivors, hazard ratio, p-value and BI
result

@param nTrials The number of virtual trials to be simulated

@param nCases The number of cases in every virtual trial

@param mediana The median time of events in the censoring dataset

@param pCure The cure-rate in the complete follow-up dataset

@param pCensoring The proportion of censored cases

@param bias The bias function

@return a data.frame

simulTrials <- function(nTrials = 300, nCases = 1000, mediana = 20, pCure = 0,

pCensoring = "runif(1, 0.05, 0.95)", ...) {

df <- data.frame(matrix(nrow = 3 * nTrials, ncol = 9), stringsAsFactors =
FALSE)

colnames(df) <-

c("type", "nCases", "pCensored", "pCured", "hr", "pValue", "SQBI", "UMBI",
"SABI")

for(i in 1:nTrials) {

nCases1 <- eval(parse(text = nCases))

median1 <- eval(parse(text = mediana))

pCure1 <- eval(parse(text = pCure))

pCensor1 <- eval(parse(text = pCensoring))

group0 <- completeFA(nCases = nCases1, median1, 100, pCure = pCure1)

timeCensored <- timeCensoring (group0, Q_A = median1, 100, pCensoring =
pCensor1)

interimCensored <- interimCensoring(group0, Q_A = median1, 100, pCensoring =
pCensor1)

caseCensored <- caseCensoring (group0, pCensoring = pCensor1)

df[3 * i - 2:0, "type"] <- c("time", "interim", "case")

df[3 * i - 2:0, -1] <- rbind(trialResults(group0, timeCensored),

trialResults(group0, interimCensored),

trialResults(group0, caseCensored))

}

df

}

This takes a cancer clinical dataset and returns the row of a data.frame with
columns: name of the trial, number of cases, proportion of censoring, SQBI and
SABI indexes and reference

@param db The clinical cancer dataset

@param trial The name of the trial

@param reference The bibliographic reference

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@return a 1 row data.frame

clinicalBias <- function(db, trial, reference = "", time = "time", status =
"status", event = 1) {

events <- db[, status] == event

censor <- db[, status] != event

lastEvent <- max(db[events, time], na.rm = FALSE)

bi1 <- SQBI(db, time = time, status = status, event = event)

bi3 <- SABI(db, time = time, status = status, event = event)

data.frame(trial = trial,

n = nrow(db),

pCens = sum(censor) / nrow(db),

SQBI = bi1,

SABI = bi3,

reference = if(reference == "") "" else paste0("\\\cite{", reference, "}"),

row.names = "", stringsAsFactors = FALSE)

}

## Plot functions

This function plots a Kaplan-Meier curve from a survival dataset. If the
logical newplot is TRUE, the curve is plotted from scratch, otherwise, the
curve is added to the plot

@param dataset The survival dataset with columns c(time, status, x)

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@param x Labels for the groups in the dataset

@param newplot A logical if the curve is to be drawn from scratch or added

@param col The color for the curve

@return a plot A Kaplan-Meier curve with the median

KMplot <- function(dataset, time = "time", status = "status", event = 1, x =
"x",

newplot = FALSE, col = "black", ...) {

op <- par(mar = c(3, 4, 4, 0.8), mgp = c(2, 0.6, 0))

if(newplot) {

xlim = c(0, max(dataset$time, na.rm = FALSE) * 1.1)

plot(NA, type="n", xlab = "Time", ylab = "Proportion", las = 1,

xaxs = "i", yaxs = "i", xlim = xlim, ylim = c(1e-2, 1), ...)

abline(h = 1:3 * 0.25, lty = "dotted")

}

formula1 <- Surv(time, status == event) ~ x

survfit1 <- survfit(formula1, data = dataset)

median1 <- summary(survfit1)$table["median"]

lines(survfit1, mark.time = TRUE, conf.int = FALSE, col = col, ...)

lines(c(median1, median1), c(-1, 0.5), lty = "dotted", lwd = 2, col = col)

par(op)

}

This function takes the complete follow-up and censored datasets and calls
KMplot to plot the Kaplan-Meier curves. The censored dataset is plotted in
black and the censored dataset in color (red for _time censoring_ , blue for
_interim censoring_ , and green for _case censoring_). In addition, the
censored dataset is split into two datasets: events dataset (plotted in
orange), and censored in the events interval (plotted in purple). The hazard
ratio and p-value comparing the two datasets are written in the title after a
label

@param group0 The complete follow-up dataset

@param group1 The censored follow-up dataset

@param time The name of the column corresponding to time

@param status The name of the column correspondint to status

@param event The label of status corresponding to censoring

@param x Labels for the groups in the dataset

@param newplot A logical if the curve is to be drawn from scratch or added

@param main Is a label to compose the plot title

@param col The color for the censored curve

@return a plot The Kaplan-Meier curves of complete follow-up and censored
datasets

trialPlot <- function(group0, group1, time = "time", status = "status", event
= 1, x = "x",

main = "", col = "black", ...) {

KMplot(group0, newplot = TRUE, lwd = 2)

KMplot(group1, col = col, ...)

lastEvent <- max(group0[group0$status == event, time])

events <- group1[group1[, status] == event, ]

censor <- group1[group1[, status] != event, ]

censor$status <- 1

KMplot(events, col = "darkorange", lwd = 2)

aCensor<- censor[censor[, time] < lastEvent, ]

KMplot(aCensor, col = "purple", lwd = 2)

censUncens <- rbind(group1, group0)

censUncens$x <- factor(censUncens$x, levels = c(group0$x[[1]], group1$x[[1]]))

cox <- lapply(getHR(censUncens), signif, digits = 3)

title(main = paste0(main, ", hr: ", cox$hr, " p: ", cox$p), ...)

subgroups <- c("uncensored dataset", "censored dataset", "events", "censored
cases")

legend("topright", bty = "n", lwd = 2, xpd=TRUE, legend = subgroups,

col = c("black", col, "darkorange", "purple") )

}

This function plots a set of 3 virtual trials from a complete follow-up
dataset and the three types of censoring with the indicated parameters

@param nCases The number of cases in every dataset

@param pCensoring The proportion of censoring in the censoring dataset

@param Q99 The quantile 99 of the events distribution

@param pCure The proportion of cured/long-term survivors

@param mediana The median of the events distribution

@return 3 plots Calls trialPlot 3 times

virtualTrials <- function(nCases = 1000, pCensoring = 0.5, Q99 = 100, pCure =
0.5, mediana = 25) {

group0 <- completeFA(nCases, mediana, Q99, pCure = pCure)

group1 <- timeCensoring(group0, mediana, Q99, pCensoring = pCensoring)

group2 <- interimCensoring(group0, mediana, Q99, pCensoring = pCensoring)

group3 <- caseCensoring(group0, pCensoring = pCensoring)

trialPlot(group0, group1, main = "Time censoring", col = "red", line = -1)

trialPlot(group0, group2, main = "Interim censoring", col = "blue", line = -1)

trialPlot(group0, group3, main = "Case censoring", col = "green", line = -1)

}

This function plots a virtual experiment as a scatter plot of _hazard ratio_
vs _proportion of censoring_ where every experiment is a dot with the
corresponding color to the type of censoring (red for _time censoring_ , blue
for _interim censoring_ , and green for _case censoring_)

@param df The experiment dataframe. Every row, a virtual trial, is represented
as a dot

@param cols The colores to be used for the 3 types of censoring

@return a scatter plot

censorPlot <- function(df, cols = c("red", "blue", "green"), ...) {

types <- c("time", "interim", "case")

pch <- ifelse(df$pValue < 0.05, 19, 21)

col <- cols[factor(df$type, levels = types)]

with(df, plot(pCensored, hr, las = 1, pch = pch, col = col, ylim = c(0, 1.2),
...))

caseCensored <- df[df[, "type"] == "case", ]

lm1 <- lm(formula = hr ~ pCensored, data = caseCensored)

abline(h = 1, lm1)

df$sig <- (df$pValue < 0.05) * 1

leg <- legend("bottomleft", bty = "n", paste(types, "censoring"))

points(x = leg$text$x + leg$rect$w, y = leg$text$y, col = cols, pch = 19)

points(x = leg$text$x + leg$rect$w * 1.3, y = leg$text$y, col = cols)

text (x = leg$text$x - .05, y = leg$text$y[1] + 0.05, "p-value <0.05 >0.05",

pos = 4, offset = 6, ...)

legend("bottomright", bty = "n",

legend = paste("r =", with(caseCensored, signif(cor(hr, pCensored), 3))))

}

This function plots a virtual experiment as a scatter plot of _hazard ratio_
vs _bias index_ where every experiment is a dot with the corresponding color
to the type of censoring (red for _time censoring_ , blue for _interim
censoring_ , and green for _case censoring_)

@param df The experiment dataframe. Every row, a virtual trial, is represented
as a dot

@param cols The colores to be used for the 3 types of censoring

@return a plot A scatter plot of _hazard ratio_ vs _proportion of censoring_.
Every virtual trial is a dot with the corresponding color to the type of
censoring (red for _time censoring_ , blue for _interim censoring_ , and green
for _case censoring_)

biasPlot <- function(df, cols = c("red", "blue", "green"), BI, ...) {

X <- df[, BI]

types <- c("time", "interim", "case")

pch <- ifelse(df$pValue < 0.05, 19, 1)

col <- cols[factor(df$type, levels = types)]

with(df, plot(X, hr, las = 1, pch = pch, col = col, xlab = BI, ylim = c(0,
1.6), ...))

leg1 <- legend("topleft", " p-value <0.05 >0.05", bty = "n")

leg <- legend(x = leg1$rect$left, y = leg1$text$y, bty = "n", paste(types,
"censoring"))

points(x = leg$text$x + leg$rect$w, y = leg$text$y, col = cols, pch = 19)

points(x = leg$text$x + leg$rect$w * 1.3, y = leg$text$y, col = cols)

df$sig <- (df$pValue < 0.05) * 1

oCutPoints <- optimal.cutpoints(

X = BI, status = "sig", data = df, tag.healthy = 0, methods = "Youden")

sCutPoints <- summary(oCutPoints)$p.table$Global

roc <- legend("topright", bty = "n", xjust = 0.5,

legend = c("Area under ROC", sCutPoints$AUC_CI))

params <- signif(sCutPoints$Youden[[1]][1:5, ], 3)

cutoff <- params[[1]]

params <- paste(c("Cutoff", "Sensit", "Specif", "posPV", "negPV"), params, sep
= ": ")

abline(h = 1, v = cutoff)

threshold <- min(BI[df$pValue < 0.05], na.rm = FALSE)

legend("bottomleft", legend = params, bty = "n")

}

## Figures

### Figure 1

Simulation of 12 virtual experiments based on 4 complete follow-up datasets
with median survival of 25 and 50 time units ($Q_{99}=100$), and cure-rate of
0 and 0.4. For every complete follow-up dataset with 1000 cases, three
censored datasets are obtained with the three different mechanisms of
censoring and proportion of censoring of 50%.

set.seed(1963)

op <- par(mfrow = c(4, 3), cex = 0.8)

virtualTrials(pCure = 0, mediana = 25)

virtualTrials(pCure = 0, mediana = 50)

virtualTrials(pCure = 0.4, mediana = 50)

virtualTrials(pCure = 0.4, mediana = 50)

par(op)

### Virtual experiments

Simulates different experiments, every one with 3000 virtual trials based on
1000 random complete follow-up datasets with a range of median distribution of
events and cure-rates

n = 1000

set.seed(1963); experiment1 <-

simulTrials(nTrials = n, mediana = "runif(1, 5, 50)", pCure = "0")

set.seed(1963); experiment2 <-

simulTrials(nTrials = n, mediana = "runif(1, 50, 95)", pCure = "0")

set.seed(1963); experiment3 <-

simulTrials(nTrials = n, mediana = "runif(1, 5, 50)", pCure = "0.5")

set.seed(1963); experiment4 <-

simulTrials(nTrials = n, mediana = "runif(1, 50, 95)", pCure = "0.5")

set.seed(1963); experiment5 <-

simulTrials(nTrials = n, mediana = "runif(1, 5, 95)", pCure = "runif(1, 0,
0.8)")

### Figure 2

Simulation of 3000 virtual trials based on 1000 random complete follow-up
datasets with medians of 5 to 95 time units and cure-rate of 0 to 0.8; time
censoring (red), interim censoring (blue) and case censoring (green) datasets
are simulated with probability of censoring of 5% to 95%.

op <- par(cex = 0.8)

censorPlot(experiment5, xlim = c(0.075, 0.925), xlab = "proportion of
censoring",

main = "median: 5-95 pCure = 0.-0.8")

par(op)

### Figure 3

Simulation of 2 experiments, each of them with 3000 trials based on 1000
complete follow-up datasets and the 3 mechanisms of censoring described, with
random proportion of censoring in the range 0.05 to 0.95.

op <- par(mfrow = c(1, 2))

biasPlot(experiment1, BI = "SQBI", main = "median: 5-50 pCure = 0", xlim =
c(0.2, 1.9))

biasPlot(experiment5, BI = "SQBI", main = "median: 5-95 pCure = 0-0.8")

par(op)

### Figure 4

Simulation of 4 virtual experiments with a 2x2 design: right skewed Weibull
distributions or not, and cure-rate of 0 or 50%.

op <- par(mfrow = c(2, 2))

biasPlot(experiment1, BI = "UMBI", main = "median: 5-50 pCure = 0.")

biasPlot(experiment2, BI = "UMBI", main = "median: 50-95 pCure = 0.")

biasPlot(experiment3, BI = "UMBI", main = "median: 5-50 pCure = 0.5")

biasPlot(experiment4, BI = "UMBI", main = "median: 50-95 pCure = 0.5")

par(op)

### Figure 5

Simulation of 2 experiments, each of them with 3000 trials based on 1000
complete follow-up datasets and the 3 mechanisms of censoring described, with
random proportion of censoring in the range 0.05 to 0.95.

op <- par(mfrow = c(1, 2))

biasPlot(experiment1, BI = "SABI", main = "median: 5-50 pCure = 0", xlim =
c(0.2, 1.9))

biasPlot(experiment5, BI = "SABI", main = "median: 5-95 pCure = 0-0.8", xlim =
c(0.2, 2.5))

par(op)

## Clinical available datasets

Available clinical datasets were retrieved from Documentation for package
$survival$ version 2.44-1.1, based on 6 published articles
(https://stat.ethz.ch/R-manual/R-devel/library/survival/html/00Index.html).
Bias indexes SQBI and SABI were applied to the datasets and to the research
arms where available.

### Acute Myelogenous Leukemia survival data

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/aml.html

trialsBI <- clinicalBias(aml, "Acute Myelogenous Leukemia survival data",
"Miller1997")

### Bladder Cancer Recurrences

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/bladder.html

data(bladder)

bladder1$time <- with(bladder1, stop - start)

bladder1$status <- (bladder1$status == 1 | bladder1$status == 2) * 1

names(bladder1)[2] <- "x"

trialsBI <- rbind(trialsBI,

clinicalBias(bladder1, "Bladder Cancer Recurrences - bladder1", "Wei1989"))

placebo <- bladder1[bladder1$x == "placebo", ]

pyridoxine <- bladder1[bladder1$x == "pyridoxine", ]

thiotepa <- bladder1[bladder1$x == "thiotepa", ]

trialsBI <- rbind(trialsBI, clinicalBias(placebo, " - placebo", ""))

trialsBI <- rbind(trialsBI, clinicalBias(pyridoxine, " - pyridoxine", ""))

trialsBI <- rbind(trialsBI, clinicalBias(thiotepa, " - thiotepa", ""))

bladder2$time <- with(bladder2, stop - start)

names(bladder2)[2] <- "x"

names(bladder2)[7] <- "status"

trialsBI <- rbind(trialsBI,

clinicalBias(bladder2, "Bladder Cancer Recurrences - bladder2", "Wei1989"))

rx1 <- bladder2[bladder2$x == 1, ]

rx2 <- bladder2[bladder2$x == 2, ]

trialsBI <- rbind(trialsBI, clinicalBias(rx1, " - rx1", ""))

trialsBI <- rbind(trialsBI, clinicalBias(rx2, " - rx2", ""))

### NCCTG Lung Cancer Data

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/lung.html

data(lung)

lung$status <- (lung$status == 2) * 1

lung$x = ""

trialsBI <- rbind(trialsBI, clinicalBias(lung, "NCCTG Lung Cancer Data",
"Loprinzi1994"))

### Chemotherapy for Stage B/C colon cancer

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/colon.html

data(colon)

colon$x <- colon$rx

trialsBI <- rbind(trialsBI,

clinicalBias(colon, "Chemotherapy for Stage B/C colon cancer", "Moertel1990"))

Obs <- colon[colon$rx == "Obs", ]

Lev <- colon[colon$rx == "Lev", ]

LFU <- colon[colon$rx == "Lev+5FU", ]

trialsBI <- rbind(trialsBI, clinicalBias(Obs, " - Observation", ""))

trialsBI <- rbind(trialsBI, clinicalBias(Lev, " - Levamisol", ""))

trialsBI <- rbind(trialsBI, clinicalBias(LFU, " - Levamisol + 5FU", ""))

### Ovarian Cancer Survival Data

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/ovarian.html

data(ovarian)

ovarian$time <- ovarian$futime

ovarian$status <- ovarian$fustat

ovarian$x <- ovarian$rx

trialsBI <- rbind(trialsBI,

clinicalBias(ovarian, "Ovarian Cancer Survival Data", "Edmonson1979"))

rx1 <- ovarian[ovarian$rx == 1, ]

rx2 <- ovarian[ovarian$rx == 2, ]

trialsBI <- rbind(trialsBI, clinicalBias(rx1, " - rx1", ""))

trialsBI <- rbind(trialsBI, clinicalBias(rx2, " - rx2", ""))

### Veterans’ Administration Lung Cancer Study

https://stat.ethz.ch/R-manual/R-devel/library/survival/html/veteran.html

data(veteran)

veteran$x <- veteran$trt

trialsBI <- rbind(trialsBI,

clinicalBias(veteran, "Veterans’ Administration Lung Cancer Study",
"Kalbfleisch1980"))

standard <- veteran[veteran$trt == 1, ]

test <- veteran[veteran$trt == 2, ]

trialsBI <- rbind(trialsBI, clinicalBias(standard, " - standard", ""))

trialsBI <- rbind(trialsBI, clinicalBias(test, " - test", ""))

### Table 1

table1 <- trialsBI

rownames(table1) <- NULL

table1[, c("pCens", "SQBI", "SABI")] <- signif(trialsBI[, c("pCens", "SQBI",
"SABI")], 2)

highlight <- function(values, threshold = 1) {

selected <- values > threshold

values[selected] <- paste0("\\\textbf{", values[selected], "}")

values

}

table1[, "SQBI"] <- highlight(table1[, "SQBI"], 1)

table1[, "SABI"] <- highlight(table1[, "SABI"], 1)

options(xtable.comment = FALSE)

xt <- xtable(table1, format = "latex", digits = 2,

caption = "Bias in available clinical datasests")

print(xt, caption.placement = "top", type = "latex", include.rownames = FALSE,

sanitize.text.function = identity)

### sessionInfo()

    ## R version 4.0.3 (2020-10-10)
    ## Platform: x86_64-apple-darwin17.0 (64-bit)
    ## Running under: macOS Catalina 10.15.7
    ##
    ## Matrix products: default
    ## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
    ## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
    ##
    ## locale:
    ## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] xtable_1.8-4           OptimalCutpoints_1.1-4 survival_3.2-7
    ## [4] knitr_1.30
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] lattice_0.20-41 digest_0.6.27   grid_4.0.3      magrittr_2.0.1
    ##  [5] evaluate_0.14   rlang_0.4.9     stringi_1.5.3   Matrix_1.2-18
    ##  [9] rmarkdown_2.5   splines_4.0.3   tools_4.0.3     stringr_1.4.0
    ## [13] xfun_0.19       yaml_2.2.1      compiler_4.0.3  htmltools_0.5.0


# DeltaZ: An Accessible Compliant Delta Robot Manipulator for Research and
Education

Sarvesh Patil∗1, Samuel C. Alvares∗2, Pragna Mannam1, Oliver Kroemer1, F.
Zeynep Temel1
∗Indicates equal contribution.1S. Patil, P. Mannam, O. Kroemer, and F. Z.
Temel are with the Robotics Institute, Carnegie Mellon University, Pittsburgh,
PA 15213, USA {sarveshp, pmannam, okroemer, ztemel} @andrew.cmu.edu2S.C.
Alvares is with the Robotics Institute Summer Scholars at Carnegie Mellon
University, Pittsburgh, PA 15213, USA and also with Rose-Hulman Institute of
Technology, Terre Haute, IN 47803, USA<EMAIL_ADDRESS>https://github.com/ZoomLabCMU/DeltaZ

###### Abstract

This paper presents the DeltaZ robot, a centimeter-scale, low-cost, delta-
style robot that allows for a broad range of capabilities and robust
functionalities. Current technologies allow DeltaZ to be 3D-printed from soft
and rigid materials so that it is easy to assemble and maintain, and lowers
the barriers to utilize. Functionality of the robot stems from its three
translational degrees of freedom and a closed form kinematic solution which
makes manipulation problems more intuitive compared to other manipulators.
Moreover, the low cost of the robot presents an opportunity to democratize
manipulators for a research setting. We also describe how the robot can be
used as a reinforcement learning benchmark. Open-source 3D-printable designs
and code are available to the public.

###### Index Terms:

Education Robotics, Parallel Robots, Compliant Joints and Mechanisms, Additive
Manufacturing, Soft Robot Applications, Flexible Robotics, Kinematics

##  I Introduction

As the manipulation capabilities of robots increase, new application domains
will continue to emerge in semi-structured environments such as homes,
warehouses, and hospitals. Given this potential impact on daily living, it is
important that the field of robotic manipulation becomes accessible to a
broader range of students and researchers. This need for accessibility is
further strengthened by manipulation challenges becoming popular test beds for
machine learning algorithms, which are already ubiquitous. A key barrier to
exploring robotic manipulation is accessibility to hardware. Low-cost
educational robots often focus on mobility rather than manipulation, and low-
cost hardware for research is often still thousands of dollars and thus poses
a significant barrier.

To democratize robotic manipulation for hands-on education and research, we
propose a low-cost, open-source, 3D-printed delta robot design, DeltaZ, shown
in Fig. 1. This robot costs around $50$ USD and its design allows for it to be
easily assembled and controlled. A core part of the design is the three 4-bar
mechanism links and articulated platform that are 3D printed as a single piece
using a soft material such as thermoplastic polyurethane (TPU). This structure
allows the robot to safely interact with its surroundings while still
achieving high repeatability. The delta kinematics also allow the robot to
directly translate the end-effector in 3D Cartesian space, with only small
rotations, which allows the three degrees-of-freedom (DoF) robot to be used
for a wide range of manipulation tasks. Being able to easily control the 3D
position also provides a more intuitive experience for novice students when
compared to three rotational DoF serial kinematics. The compliant delta
structure also allows the robot to be used for dynamic tasks, such as striking
a ball or puck. The base plate can be easily replaced with different designs
to create various task-specific environments and structures, and additional
sensors can be easily incorporated to augment the robot’s capabilities.

Figure 1: The DeltaZ and some of its functionalities. DeltaZ is a low-cost
manipulation robot with 3-D printed mechanical parts and a compliant end-
effector (a). DeltaZ’s hook end-effector picking up a weight (b), with a stamp
tool and an ink pad, stamping a letter (c), and with a 3D-printed boot
attachment, kicks a soccer ball into the goal (d).

In the following sections, we present the design of the low-cost delta robot
and its fabrication. We evaluate the robot’s capabilities in terms of speed
and repeatability. Our experiments demonstrate the robustness of the robot and
show how it can be used as a benchmarking tool for robot reinforcement
learning, as well as a basis for multi-robot transfer learning.

##  II Related Work

###  II-A Delta Robots

Delta robots are parallel robots with 3 translational degrees of freedom where
the end-effector stays parallel to its base (as shown in Fig.3(a)). The motors
are stationed in the body which allows the light end-effectors to perform pick
and place tasks with high accuracy and precision [1]. With closed-form inverse
kinematic solutions, the precision of delta robots can be controlled
sufficiently for biomedical and surgery applications [2]. Previous work shows
3D-printed compliant end-effectors with linearly actuated delta robots can be
used for dexterous manipulation [3]. We build upon this work to make a
compliant delta robot manipulator with revolute actuation for safe
interactions with objects and people.

###  II-B Educational Robots

A number of educational robotics kits have been developed in recent years [4].
These robots are often designed to teach students about coding and make the
sense-plan-act loop accessible and engaging to K-12 students. However, these
robots tend to emphasize mobility rather than manipulation, with manipulation
often being a 1 DoF extension to the robot like a fork lift. The focus on
mobility also means that the lesson is often on the robots avoiding contacts,
or only allowing for a narrow set of interactions, rather than using
compliance to exploit contacts with the environment. Our DeltaZ robot was
designed to afford safe compliant interactions. Future goals include
developing the teaching material that one would need for a full educational
kit and course.

###  II-C Benchmarking for Robot Manipulation

Developing benchmarking tasks for robot manipulation research is challenging,
and a number of workshops at international conferences have been dedicated to
this issue [5, 6, 7]. A core problem is the need to have similar objects and
hardware to compare different algorithms without having to repeat entire
experiments. Object sets and shared experimental protocols have helped to make
experiments more reproducible [8]. However, differences in robot hardware and
low-level control can still have a significant effect on results.

To avoid hardware issues, simulators have become ubiquitous tools for
evaluating robot manipulation algorithms [9]. However, modern simulation
engines still exhibit a large simulation-to-reality gap and often provide an
over-idealized version of tasks. The simulation-to-reality gap is exacerbated
when modeling the complexities of manipulating non-rigid objects. The DeltaZ
robot makes real-robot evaluations more accessible, allowing for easier hands-
on experimentation.

Remote experimentation sites, wherein the robots are controlled remotely over
the internet, provide access to robots for a larger population of researchers
[10]. These sites allow state of the art robots to be accessed from around the
globe. However, the remote nature of these tasks makes it difficult to create
new task environments. For the DeltaZ robot, new task environments and objects
can be 3D printed given its scale. The DeltaZ also encourages hands-on
research which can often lead to new insights when observing experiments in
person.

It should also be noted that simulations and remote experimentation sites are
favored by some researchers, and especially non-robotics researchers, as they
avoid the need for setting up and maintaining a physical robot system. These
processes can be time consuming and often require additional engineering
expertise that may not be present in some labs. The DeltaZ robot provides a
simple tabletop environment that is largely self contained and easy to
assemble and maintain.

###  II-D Low-cost Research Manipulators

The need for accessible research robots has resulted in a number of low-cost
robots being developed in recent years. For example, the Locobot incorporates
a 5 DoF robot arm, and the Dynamixel Claw includes finger-like 3 DoF
manipulators[11, 12]. These manipulators use serial kinematic designs, like
traditional robot arms. Serial kinematics result in each motor affecting the
end-effector’s position and orientation. A robot will therefore often need to
have additional degrees of freedom simply to maintain a certain orientation,
or the task will need to be designed to reduce the effects of the rotations
(e.g., using a spherical end-effector). The delta design allows us to create
and perform a wide range of translational manipulation tasks with a simple 3
DoF design.

A serial design also means that the servos need to be strong enough to support
the other motors in the chain. This requirement does not only increase the
cost of the motors, but can also be taxing on the motors over time. The
parallel design of the delta robots allows us to use lower-cost servo motors.

The cost of typical robotic manipulators is in the range of thousands of
dollars [13, 14, 15]. To make the DeltaZ accessible to students and
researchers, we designed the 3D-printed revolute delta robot to cost around
$50$ USD.

##  III DeltaZ Robot Design

The core design goals of the DeltaZ are to make a robot that is precise,
versatile, low-cost, and can withstand impacts and obstacles due to its
compliance. The affordable robot enables a multitude of manipulation tasks to
be achieved and is ideal for an educational or research setting. Mechanical
design, electronics and control, as well as examples of robot functionalities
are discussed in this section.

###  III-A 3D-Printing and Component Overview

DeltaZ is made of a total of 42 pieces, including individual screws. An
overview of the individual pieces is shown in Fig.2, and assembled and labeled
in Fig.3. The black pieces (components 1-5 and 14 in Fig.2) are all 3D printed
from Polylactic Acid (PLA). These components include a housing to encase the
electronics and mount the servos, as well as a set of legs and a base plate
for supporting the delta robot’s body when used in a top-down manner.

The mechanical assembly using screws and bolts provides modularity, i.e., all
of the components can be replaced if needed as nothing is glued. The time to
assemble the robot varies based upon the users experience level, but can
generally be done within an hour. A video tutorial on DeltaZ’s website [16]
walks users through the assembly process.

All designs are open hardware which allows communities of users to share ideas
and collaborate on adapting the overall design and functionality of the robot
to specific applications [17].

Figure 2: All of the parts required to build the robot, including the
3D-printed parts. Parts include (1) Body, (2) Base, (3) Cap, (4) Leg, (5)
Divider, (6) Expansion Board, (7) Nano, (8) 8 mm Screw, (9) 10 mm Screw, (11)
Nut, (12) Self Tapping Screw, (13) Servo Screw, (14) Forearm, (15) Servo, (16)
Compliant End-Effector. Note that Self Tapping Screws (12) and the Servo
Screws (13) are included with the Servos (15). Figure 3: (a) Assembled DeltaZ.
The major parts include the body, servo motor, forearm, compliant end-effector
(printed with white material), legs, and base. The body encloses all motors
and electronics and is supported by three legs which form a tripod. (b) The
compliant end-effector of DeltaZ can be 3D-printed with two orthogonal
revolute joints printed as living hinges. Parameters like k and L can be
modified by users. The four central holes to attach a variety of end-effectors
using M3 screws. This component bolts to the rigid forearms via the outer
three screw holes. (c) DeltaZ’s allowable workspace shown as a transparent
cylindrical solid with a diameter of 60mm and height of 75 mm.

###  III-B Soft Mechanisms and Living Hinges

The white component (16) in Figs. 2, 3(a), and 3(b) is referred to as the
compliant end-effector, as it is printed from a flexible material, such as TPU
95A. Delta robots use a set of parallel 4-bar mechanisms (parallelograms) to
maintain the end-effector’s orientation. This part of the design often results
in a large number of additional components that increase the overall
complexity of the robot’s design. To provide a design that is simple, yet
precise, we 3D print the entire structure as a single print.

The flexible joints are made as living hinges to mimic a single revolute
joint. A combination of two living hinges performs similarly to a universal
joint, which is typically present in a conventional delta robots. Living
hinges are made by locally reducing material thickness to a compliant
articulated joint. For experiments, we printed delta parts using Ultimaker TPU
95A, PP (polypropelene), and PLA (polylactic acid) material. From most to
least rigid, the tensile moduli are 2,346.5 MPa (using ISO 527), 220 MPa
(using ISO 527), 26 MPa (using ASTM D638) for PLA, PP, and TPU, respectively
[18]. The materials were chosen for their ability to be 3D-printed and create
living hinges at the desired locations. Non-rigid materials like PP and TPU,
in addition to relatively low-torque motors, make the robot safe to users.

Design parameters like parallelogram beam and hinge thickness were chosen
based on prior work [19]. The beam thickness is $4.5$ mm for TPU and PLA
deltas and $2.5$ mm for PP delta parts as in [3, 19]. The hinges are all
$0.41$ mm thick which was as close to the desired $0.375$ mm found in previous
work but could also be printed on various 3D-printers (limited by size of
material extruded). Parameters L and k in Fig. Fig.3(b) are parallelogram link
lengths of $L=37$ mm, and $k=5.25$ mm which is the offset between the two
orthogonal revolute joints made from living hinges. These two revolute joints
should be as close as possible ($k\rightarrow 0$ mm) to minimize positional
error [3]. Using key physical parameters from previous work, we improved upon
the mechanical design of a compliant, delta-style robot by making it more
conducive for a learning and research environment. In particular, mechanical
connections were made more durable and reliable. Also, simple assembly and
fabrication processes were formalized and made open-source. Users may
experiment with the hinge thickness, parallelogram link length, and offset
between revolute joints to test the affect on kinematic behavior of the
compliant end-effector. Additionally, changes in the resulting workspace could
be helpful for specific manipulation tasks.

###  III-C Base Plate

The base plate (part 2 in Fig. 2) serves as a support surface for the robot’s
manipulation tasks. The basic plate is flat and simply ensures that the robot
cannot push off of the underlying table. However, the plate can also be easily
replaced with other base plates to create different task-specific
environments. For our experiments, one of the base plates includes a
potentiometer mounting for a dial turning experiment. Other designs could
including fixtures for different tasks. In this manner the robot can be easily
adapted to explore different tasks and provide a self-contained environment
for benchmarking tasks. Switching out the base plate can be done in a matter
of minutes by simply unscrewing the legs, exchanging the plates, and
reattaching the legs. The base plate can also be removed to allow the robot to
perform simple locomotion tasks with a single articulated leg.

###  III-D Arduino, Servos, and Sensorization

DeltaZ is driven by an Arduino Nano microcontroller and powered via USB. Delta
robots have base mounted motors and parallel geometry that allow for fast and
accurate motions with relatively small and low-cost motors [20, 2]. Thus, we
are able to articulate DeltaZ using affordable, 9-gram, metal geared, micro
servo motors. DeltaZ can be positioned by both forward and inverse kinematics
in open-loop control. The Arduino also allows for simple sensors, such as
buttons, light sensors, or potentiometers, to be easily incorporated into the
platform.

###  III-E Serial Interface

The robot is controlled externally through the serial port interface of the
Arduino. The interface allows the desired angles to be directly specified or
the desired x-y-z location of the end-effector to be given. For the latter,
the inverse kinematics are computed directly on the Arduino to compute the
corresponding desired angles based on a rigid-link model of the robot. We
limit the workspace to a cylindrical region of height $40$ mm and diameter 60
mm, as shown in Fig. 3(c), such that the end-effector cannot collide with the
legs of the robot. The interface can also be used to read the values of
sensors connected to the Arduino. All software and hardware designs are open-
source.

##  IV DeltaZ for RL Benchmarking

To show the efficacy of the DeltaZ for real-robot benchmarking, we have the
robot apply reinforcement learning to acquire a dial turning skill. The
learning processes is repeated multiple times across three different copies of
the robot to demonstrate the similar outcomes.

###  IV-A Dial Turning Task

The example benchmarking task is designed around a potentiometer mounted in
the base of the robot, as shown in Fig. 4. The goal of the task is for the
robot to use its end-effector to turn the potentiometer to match a desired
resistance value. A small lever has been attached to the potentiometer for the
robot to push against. The potentiometer is mounted in a 3D-printed base plate
that was designed for this task, and its pins are connected to the Arduino
such that its resistance values can be easily measured by the robot. This task
requires the robot to operate an articulated object through contact-based
interactions, with different amounts of force required depending on where on
the lever the robot pushes. The task was inspired by similar tasks used for
the Dynamixel claw [21]. For each epsiode, the robot receives a reward of
$R=100$ if the final angle $\phi$ is within $15^{\circ}$ of the desired angle
$\phi_{d}$, indicating a successful task completion, as well as a quadratic
cost based on the difference between the final angle and the desired angle
$R=100[\|\phi-\phi_{d}\|<15]-10^{-5}(\phi-\phi_{d})^{2}$.

Figure 4: DeltaZ robot with a potentiometer mounted on the base (a). The
potentiometer (b) is connected to an analog input pin on the Arduino.
Corresponding to the resistance value of the potentiometer, we obtain values
from 0 to 1023 as analog input from the Arduino.

To minimize the amount of human effort in running the evaluations, we
incorporate an automatic resetting procedure for re-positioning the dial
between trials. The end-effector goes behind the lever and rotates it back to
the starting angle, at which point a new episode can be executed.

The experiments were conducted across three different robots, with seven full
learning processes run on each robot. The resulting 21 trails were performed
to demonstrate the reproducibility of the learning process across different
robots.

###  IV-B Skill Parameterization

The end-effector of the DeltaZ robot is restricted to a $30$ mm radius around
the origin. Thus, in order to efficiently represent the action space, we
parameterize it in the polar $(\rho,\theta)$ domain, where $\rho$ is the
radius and $\theta$ is the angle from the positive x-axis. The z position of
the end-effector is fixed during the skill execution at a height that it can
push the potentiometer’s lever.

Each skill execution is then defined by two points $(\rho_{1},\theta_{1})$,
and $(\rho_{2},\theta_{2})$ that define two waypoints. The robot moves down at
the first waypoint and then moves across to the second waypoint in a straight
line. Thus, the waypoints are defined in polar coordinates, but the trajectory
itself is in a straight Cartesian x-y line.

The goal of the reinforcement learning is to learn a set of skill parameters
$(\rho_{1},\theta_{1},\rho_{2},\theta_{2})$ for achieving a desired
potentiometer angle.

###  IV-C Skill Learning

To demonstrate skill learning on the DeltaZ platform, we learn the dial
turning task using episodic Relative Entropy Policy Search (eREPS) [22]. eREPS
models the distribution over the skill parameters as a 4D Gaussian
distribution. We normalize the parameter values to be within a range of -1 to
1, and initialize the Gaussian distribution with mean $0.4$ and a diagonal
covariance matrix with non-zero elements of $0.15$.

As a model-free policy search approach, eREPS iterates between evaluating
batches of sampled parameters on the real robot and updating the Gaussian
policy based on the resulting rewards. For our experiments, the robot rolls
out 20 episodes initially, and then 10 episodes for subsequent iterations
between each policy update. The actions sampled from the Gaussian distribution
are passed to the Arduino which moves the end-effector to the desired
location, reads the resistance value of the potentiometer, and sends it back
to the computer to compute the reward. Given the automatic resets and
potentiometer recording, the data collection process can be run fully
autonomously.

The policy updates of eREPS attempt to maximize the expected return while
limiting the KL divergence between the old policy and the new policy. For our
updates, the robot utilizes the most recent 20 samples for each policy update,
i.e., from the last two policies. Retaining samples from multiple previous
policies is a common practice for eREPS to further improve learning stability.
The eREPS algorithm is terminated when all 10 trials in the previous batch
were successful at completing the task.

Figure 5: Convergence of parameters for a single Robot. Figure 6: Gaussian
plots of rewards across different robots.

##  V Evaluations

Our experiments evaluated the performance of the robot performing tablet
interaction tasks and the dial-turning RL benchmarking task. The results of
our evaluations are given in this section.

###  V-A Drawing Task Evaluations

We used the DeltaZ robot to perform a series of drawing tasks to characterize
its workspace. Two end-effectors of TPU and PP are used for the experiments.
The end-effectors were equipped with a capacitive stylus pen and the drawing
was performed on a tablet, as shown in Fig. 7 (a). To test the repeatability
of the DeltaZ, lines along x and y axes were drawn at different lateral and
vertical distances. First, the robot is commanded to make a plus sign at the
origin of the plane. Then, the stylus pen followed a straight trajectory along
x-axis with three different y-values (Fig. 7 (b)). The same experiments are
repeated to observe the accuracy and precision along the y-axis (Fig. 7 (c)).
The results over 10 repetitions show that the accuracy of the robot is not
high. However, it performs the tasks with very high precision.

The comparison between end-effectors printed using different materials is
shown in Fig. 8 following a circular trajectory with different radii and
different vertical distances. Both versions of DeltaZ were able to follow a
fairly precise trajectory with a changing radius. However, as we increased the
vertical distance and pushed the stylus pen into the tablet, we observed that
the trajectory got distorted as the vertical distance of the TPU end-effector
decreased. However, DeltaZ with the PP end-effector was not able to finish the
tasks at different values along z-axis, due to the lack of compliance.

In order to measure the maximum speed of DeltaZ, we performed experiments
where the end-effector follows a straight line of 40mm along the x-axis
multiple times. The maximum speed at which DeltaZ was still able to follow the
full trajectory is calculated as 0.17 m/s.

Figure 7: DeltaZ is placed on an tablet with a stylus pen to characterize the
workspace of the robot (a). The accuracy and precision of DeltaZ over 10
trials for straight lines along x (b) and y (c) axes are shown.  Figure 8: The
comparison between end-effectors that are printed with TPU (a and b) and PP (c
and d). Radius of the circle varies between 30 mm and 10 mm, colors indicating
different radii (a and c). As the distance between the stylus pen and the
touchpad reduces, the shape changes due to the compliance of the end-effector
(b and d). The effect of the material is clearly observed. Since the PP end-
effector cannot conform to the environment as well as TPU end-effector, it
snapped off of the forearms, hence resulting with a non-uniform shape (green
line).

###  V-B RL Benchmarking

The eREPS algorithm was run on three different robots seven times each for
learning the dial turning task. The goal of this experiment is to evaluate the
reproducability of the results across the different robots in order to
demonstrate their utility for benchmarking.

The results of the experiments are shown in Fig. 5 and 6. Fig. 5 shows an
example of the average rewards and the distribution over the skill parameters
for each policy update of one of the 21 learning processes. Fig. 6 shows the
distribution over the seven average rewards for each of the robots. The error
bars correspond to two standard errors.

As shown in Fig. 5, the robot starts with a broad Gaussian distribution to
explore the parameter space. This initial exploration resulted in the end-
effector colliding against the potentiometer axle several times. Due to the
compliance of the end-effector, these collisions did not cause any damage and
the robot could continue to perform the entire training process without human
intervention.

The main result of this experiment is the similarity of the learning curves
across the different robots. The means are close together, with a large amount
of overlap between the standard error regions. This indicates that the
evaluations performed on the different robots are comparable. These results
demonstrate that the three different robots could be used for benchmarking
algorithms. The ability to reproduce the task environment and run the
experiment autonomously also allows for easier reproduction of results across
robots.

As an additional evaluation, the learned skills were executed on the other
robots, and we found that each of them succeeded using direct zero-shot
transfer. This result demonstrates the potential of using the DeltaZ robots
for multi-robot training and robot-robot transfer research in the future.

##  VI Conclusion

We proposed the DeltaZ design for creating an accessible robot platform for
research and education. The 3D printed robot is easy to assemble and maintain,
and it can be used to perform a variety of translational manipulation tasks.
Our experiments demonstrated the robot’s ability to perform precise repetitive
tasks that exploit the robot’s inherent compliance. We also demonstrated the
utility of the robot as a benchmarking platform. In the future, we will
explore extending the platform with additional sensors to perform various
tasks, as well as create accompanying teaching materials to develop an
educational kit.

##  VII Acknowledgements

The research presented in this work was funded by the National Science
Foundation under project Grant No. CMMI-2024794.

## References

  * [1] J.-P. Merlet, _Parallel robots_. Springer Science & Business Media, 2005, vol. 128.
  * [2] H. McClintock, F. Z. Temel, N. Doshi, J.-s. Koh, and R. J. Wood, “The milliDelta: A high-bandwidth, high-precision, millimeter-scale Delta robot,” _Science Robotics_ , vol. 3, no. 14, Jan. 2018, publisher: Science Robotics Section: Research Article. [Online]. Available: https://robotics.sciencemag.org/content/3/14/eaar3018
  * [3] P. Mannam, A. Rudich, K. Zhang, M. Veloso, O. Kroemer, and F. Temel, “A low-cost compliant gripper using cooperative mini-delta robots for dexterous manipulation,” _Robotics: Science and Systems XVII_ , 2021.
  * [4] T. Sapounidis and D. Alimisis, “Educational robotics for STEM: A review of technologies and some educational considerations,” Sept. 2020, pp. 167–190.
  * [5] B. Calli, A. Walsman, A. Singh, S. Srinivasa, P. Abbeel, and A. M. Dollar, “Benchmarking in manipulation research: Using the yale-cmu-berkeley object and model set,” _IEEE Robotics & Automation Magazine_, vol. 22, no. 3, pp. 36–52, 2015.
  * [6] “Cfp: Benchmarking protocols for robotic manipulation.” [Online]. Available: https://www.ieee-ras.org/publications/ra-l/special-issues/benchmarking-protocols-for-robotic-manipulation
  * [7] “Iros 2021 workshop - benchmarking of robotic grasping and manipulation: protocols, metrics and data analysis.” [Online]. Available: https://www.robotics.qmul.ac.uk/events/iros-2021-workshop/
  * [8] B. Calli, A. Walsman, A. Singh, S. Srinivasa, P. Abbeel, and A. M. Dollar, “Benchmarking in manipulation research: Using the yale-cmu-berkeley object and model set,” _IEEE Robotics Automation Magazine_ , vol. 22, no. 3, pp. 36–52, 2015.
  * [9] T. Yu, D. Quillen, Z. He, R. Julian, K. Hausman, C. Finn, and S. Levine, “Meta-world: A benchmark and evaluation for multi-task and meta reinforcement learning,” _CoRR_ , vol. abs/1910.10897, 2019. [Online]. Available: http://arxiv.org/abs/1910.10897
  * [10] M. Chevalier-Boisvert, F. Golemo, Y. Cao, B. Mehta, and L. Paull, “Duckietown environments for openai gym,” https://github.com/duckietown/gym-duckietown, 2018.
  * [11] A. Murali, T. Chen, K. V. Alwala, D. Gandhi, L. Pinto, S. Gupta, and A. Gupta, “Pyrobot: An open-source robotics framework for research and benchmarking,” _arXiv preprint arXiv:1906.08236_ , 2019.
  * [12] “Dexterous manipulation with reinforcement learning: Efficient, general, and low-cost.” [Online]. Available: https://bair.berkeley.edu/blog/2018/08/31/dexterous-manip/
  * [13] H. Zhu, A. Gupta, A. Rajeswaran, S. Levine, and V. Kumar, “Dexterous manipulation with deep reinforcement learning: Efficient, general, and low-cost,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 3651–3657.
  * [14] “Thinkbot solutions: Onrobot.” [Online]. Available: https://thinkbotsolutions.com/collections/onrobot
  * [15] “Thinkbot solutions: Robotiq.” [Online]. Available: https://thinkbotsolutions.com/collections/robotiq
  * [16] Sam Alvares, “Building DeltaZ,” Aug. 2021. [Online]. Available: https://www.youtube.com/watch?v=Mg2i0toDOXg
  * [17] “Open Source Hardware Association -.” [Online]. Available: https://www.oshwa.org/
  * [18] “Ultimaker support: Materials.” [Online]. Available: https://support.ultimaker.com/hc/en-us/categories/360002336619-Materials
  * [19] P. Mannam, O. Kroemer, and F. Z. Temel, “Characterization of Compliant Parallelogram Links for 3D-Printed Delta Manipulators,” in _Experimental Robotics_ , B. Siciliano, C. Laschi, and O. Khatib, Eds. Cham: Springer International Publishing, 2021, vol. 19, pp. 75–84, series Title: Springer Proceedings in Advanced Robotics. [Online]. Available: https://link.springer.com/10.1007/978-3-030-71151-1_7
  * [20] M. López, E. Castillo, G. García, and A. Bashir, “Delta robot: Inverse, direct, and intermediate Jacobians,” _Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science_ , vol. 220, no. 1, pp. 103–109, Jan. 2006. [Online]. Available: http://journals.sagepub.com/doi/10.1243/095440606X78263
  * [21] M. Ahn, H. Zhu, K. Hartikainen, H. Ponte, A. Gupta, S. Levine, and V. Kumar, “ROBEL: RObotics BEnchmarks for Learning with low-cost robots,” in _Conference on Robot Learning (CoRL)_ , 2019.
  * [22] J. Peters, K. Mülling, and Y. Altün, “Relative entropy policy search,” in _Proceedings of the Twenty-Fourth AAAI Conference on Artificial Intelligence_ , ser. AAAI’10. AAAI Press, 2010, p. 1607–1612.
  * [23] “Craftcloud® by All3DP | Your 3D Printing Service Marketplace.” [Online]. Available: https://craftcloud3d.com/