R: Clinical Trial Simulator

The first version of the clinical trial simulator is online now.

It uses a dose-response model based on the Hill equation which basically is a saturation model. The assumption is that a cell is covered by a finite number of receptors. When the treatment is applied, a fraction of the receptors (depending on the dosage) is interacting with the treatment molecule and every interaction leads to a response. When every receptor is binding to the treatment, the response saturates.

It's based on the MSToolkit package and is hosted online at shiny apps under

https://mzhku.shinyapps.io/trialsim/

R: Sample Size Calculator

Given a beta, how many patients need to be enrolled in a two-armed study of conventional against experimental therapy? Significance is set to 0.05, recruitment period and follow-up period can be defined together with median disease free survival times for the two study groups (implemented in R/Shiny):

Github repository

Try it (also mobile)!
https://mzhku.shinyapps.io/size/

R Shiny: Some reactivity by example

Start R from the directory containing ui.r and server.r, import the shiny library and start the (development) server using runApp().

server.r

function(input, output) {

    rv <- reactiveValues()

    # Only modify and output reactive input

    output$numOut1 <- renderText({

        input$n+10

    })

    

    # Modify, reassign and output reactive input.

    output$numOut2 <- renderText({

        rv$a <- input$n+20

    })

    

    # Use modified reactive input, modify and output.

    output$num1 <- renderText({

        rv$s <- rv$a+runif(1)

    })

    

    # Assign reactive input locally, modify locally, print.

    output$num2 <- renderText({

        k <- input$n

        paste(as.character(runif(1)), "here", k-10)

    })

    

    # Use reactive value from another reactive function,

    # modify and print.

    output$let <- renderText({

        paste(as.character(rv$s), sample(letters, 1), sep="")

    })

}

ui.r

fluidPage(

    numericInput(inputId="n", "Shiny Reactivity Example", value=25),

    textOutput(outputId="numOut1"),

    

    textOutput(outputId="numOut2"),

    

    textOutput(outputId="num1"),

    textOutput(outputId="num2"),

    textOutput(outputId="let")

)   

R: Plotting multiple categorical data

Frequently, values are plotted against a categorical argument. The categorical argument can be thought of something like a weekday label or a step in a specific technical procedure ("Step 1", "Step 2", ...). Of course, these arguments could be converted to a numeric value used for plotting and after plotting, one would overwrite the x-axis labels to get the categorical data on the axis. But if Excel can do it, then it should be possible to do it in R as well. Here's how.

The data in the example is from the book "Design and Analysis of Clinical Trials: Concepts and Methodologies" by Chow (full data set).

Here, the categorical data (x-axis) are five repeated measurements of diastolic blood pressure, denoted "DBP1", ..., "DBP5", where DBP1 is the baseline value. 40 subjects have been treated with either a hypertensive agent (A) or placebo (B). After averaging over treatment groups we end up with two groups we want to inspect.

> d <- aggregate(dat[, 3:7], by=list(TRT=dat$TRT), FUN=mean)

> d
  TRT   DBP1  DBP2   DBP3   DBP4   DBP5
1   A 116.55 113.5 110.70 106.25 101.35
2   B 116.75 115.2 114.05 112.45 111.95

How can now d be plotted to immediately make the difference between the treatments visible? Lattice can plot multiple data records (the A- and B-rows) against categorical data, but for that, the data has to be in long format. This can be generated using the melt function from the reshape2 package.

> m <- melt(d, id="TRT", measure.vars=colnames(d[, 2:ncol(d)]))

> m
   TRT variable  value
1    A     DBP1 116.55
2    B     DBP1 116.75
3    A     DBP2 113.50
4    B     DBP2 115.20
5    A     DBP3 110.70
6    B     DBP3 114.05
7    A     DBP4 106.25
8    B     DBP4 112.45
9    A     DBP5 101.35
10   B     DBP5 111.95

The melted data.frame m can now be plotted.

> xyplot(m$value~m$variable, type="o", group=m$TRT, auto.key=list(T))

The resulting figure looks like this:

 

Clustering data with k-means and plotting for exploratory analysis

I clustered data from my paper A computational method for the systematic screening of reaction barriers in enzymes: searching for Bacillus circulans xylanase mutants with greater activity towards a synthetic substrate using k-means clustering and plotted the results using the `Lattice` library in R.

This analysis could be used to further automate selection of reaction barriers for including or excluding from analysis if the barrier profile is not physically meaningfull.

Raw data:
http://pastebin.com/41TrNNvH

A detailed report on the analyis is available here:
http://rpubs.com/mzh/167769

Github repository:
https://github.com/mzhKU/Enzyme-Barrier-Clustering

The critical part of the code is where the two plots (the k-means clusters, red lines, and the raw data of each cluster, blue lines) are overlayed using the lattice panel plot construct (in 003_execute.r):

# 'group' required to prevent drawing of a continuous line connecting the last
# data point of a mutant with the first data point of the next mutant barrier.
xyplot(Barrier_Cl ~ x|Cluster_fac, data=kmdf, iedf_i=iedf,
       ylim=c(-20, 50), ylab="Reaction Energy [kcal/mol]",
       xlab="Reaction Coordinate", xlim=c(1, 12), strip=T,
       # Note:
       # 'x' and 'y' are the data from the cluster data frame 'kmdf'.
       # 'iedf_i' is the 'iedf'-data frame provided to the 'i'-nner panel function.
       panel = function(x, y, iedf_i)
       {
            panel.xyplot(x=iedf_i[iedf_i$Cluster==panel.number(), ]$x_axis,
                         y=iedf_i[iedf_i$Cluster==panel.number(), ]$Barrier_i,
                         group=iedf_i$id, subscripts=TRUE, type="o", col="blue")
            panel.xyplot(x, y, type="l", col="red", lwd=2)
       }
)

Excel like PivotTable and more with R

A table contains values as below (download from link under figure):

Item	Category	Price	Profit	Actual Profit	Calories
Beer	Beverages	$4.00	50%	$2.00	200
Soda	Beverages	$2.50	80%	$2.00	120
Chocolate Bar	Candy	$2.00	75%	$1.50	255
Ice Cream Sandwich	Frozen Treats	$3.00	67%	$2.00	240
Bottled Water	Beverages	$3.00	83%	$2.50	0
Gummy Bears	Candy	$2.00	50%	$1.00	300
Soda	Beverages	$2.50	80%	$2.00	120
Hamburger	Hot Food	$3.00	67%	$2.00	320
Popcorn	Hot Food	$5.00	80%	$4.00	500
Licorice Rope	Candy	$2.00	50%	$1.00	280
Hot Dog	Hot Food	$1.50	67%	$1.00	265
Licorice Rope	Candy	$2.00	50%	$1.00	280
Popcorn	Hot Food	$5.00	80%	$4.00	500
Popcorn	Hot Food	$5.00	80%	$4.00	500

It is very easy to summarize the data using the Excel PivotTable generator (even for multiple variables like "Profit" and "Actual Profit" simultaneously):

Example taken from the Data Smart book.

It probably doesn't make a lot of sense to summarize the data by percentage, but this is just for illustration purposes.







A similar operation is the aggregate function in R (the Concessions.xlsx data frame is loaded as con):
> aggregate(cbind(Profit, Actual.Profit) ~ Category, data=con, FUN=sum)
       Category   Profit Actual.Profit
1     Beverages 31.23333          98.5
2         Candy 23.25000          46.5
3 Frozen Treats 23.00000          69.0
4      Hot Food 45.21667         142.0



The same is possible using the xtabs function:
> xtabs(cbind(Profit, Actual.Profit) ~ Category, data=con)
              
Category           Profit Actual.Profit
  Beverages      31.23333      98.50000
  Candy          23.25000      46.50000
  Frozen Treats  23.00000      69.00000
  Hot Food       45.21667     142.00000

Note, no special operation FUN is specified in xtabs.

Occurence of items can be done when specifing "Count" in the PivotTable builder:

R can do this with the plyr package and its count function.

> count(con$Item)
                       x freq
1                   Beer   20
2          Bottled Water   13
3          Chocolate Bar   13
4  Chocolate Dipped Cone   11
5            Gummy Bears   14
6              Hamburger   16
7                Hot Dog   15
8     Ice Cream Sandwich   10
9          Licorice Rope   13
10                Nachos   15
11                 Pizza   17
12               Popcorn   16
13              Popsicle   13
14                  Soda   13




This can also be done using the extremely powerful built-in functions of the data.table type:

> con[, .N, by=Item]
                     Item  N
 1:                  Beer 20
 2:         Bottled Water 13
 3:                  Soda 13
 4:         Chocolate Bar 13
 5:           Gummy Bears 14
 6:         Licorice Rope 13
 7:              Popsicle 13
 8:    Ice Cream Sandwich 10
 9: Chocolate Dipped Cone 11
10:               Popcorn 16
11:             Hamburger 16
12:                Nachos 15
13:                 Pizza 17
14:               Hot Dog 15

Here, .N is a built-in count function applied to the data.table object.

Or with a call to table where one can break down the items by profit:

> table(con$Item, con$Profit)
                      
                        0.25 0.5 0.67 0.75 0.8 0.83
  Beer                     0  20    0    0   0    0
  Bottled Water            0   0    0    0   0   13
  Chocolate Bar            0   0    0   13   0    0
  Chocolate Dipped Cone    0  11    0    0   0    0
  Gummy Bears              0  14    0    0   0    0
  Hamburger                0   0   16    0   0    0
  Hot Dog                  0   0   15    0   0    0
  Ice Cream Sandwich       0   0   10    0   0    0
  Licorice Rope            0  13    0    0   0    0
  Nachos                   0  15    0    0   0    0
  Pizza                   17   0    0    0   0    0
  Popcorn                  0   0    0    0  16    0
  Popsicle                 0   0    0    0   0   13
  Soda                     0   0    0    0  13    0

I.e. 13 items of "Bottled Water" giving a profit of 83% each have been sold.

Excel can then apply eg. summation when breaking down items by category:

The operation to carry out (sum, count, ...) is defined when clicking on the small "i" symbol in the Values field. So in total, beer earned 80 $.

Another plyr built-in function of the data.table object allows to quickly obtain the summarized prices of the category sales, broken down by category.

> con[, sum(Price), by=list(Category, Item)]
         Category                  Item   V1
 1:     Beverages                  Beer 80.0
 2:     Beverages         Bottled Water 39.0
 3:     Beverages                  Soda 32.5
 4:         Candy         Chocolate Bar 26.0
 5:         Candy           Gummy Bears 28.0
 6:         Candy         Licorice Rope 26.0
 7: Frozen Treats              Popsicle 39.0
 8: Frozen Treats    Ice Cream Sandwich 30.0
 9: Frozen Treats Chocolate Dipped Cone 33.0
10:      Hot Food               Popcorn 80.0
11:      Hot Food             Hamburger 48.0
12:      Hot Food                Nachos 45.0
13:      Hot Food                 Pizza 34.0
14:      Hot Food               Hot Dog 22.5