A bit of background: I am a virologist and only dabble in r, so wanted to hone my skills by making a r TCID50 (tissue-culture infectious dose 50% end-point) calculator based off the Reed-Muench method (link may be paywalled). Basically what this does is calculate a proportional distance off a scored (positive/negative) matrix, and finds a point at which 50% of the values would be positive. The data looks a bit like this (note the 1/0 transitions in col 3/4):
| 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 0 | 0 | 0 | 0 |
| 1 | 1 | 1 | 1 | 0 | 0 | 0 |
| 1 | 1 | 1 | 0 | 0 | 0 | 0 |
| 1 | 1 | 0 | 0 | 0 | 0 | 0 |
| 1 | 1 | 1 | 0 | 0 | 0 | 0 |
Where each column is a step in a dilution series of the infectious agent and the 1/0 are positive/negative. I have an Excel sheet to do this too, but I like playing with r.
Actual data and packages:
library(tidyverse)
library(FSA)
Scored_Data <- structure(list(X1 = c(0, 0, 0, 0, 0, 0, 0, 0), X2 = c(0, 0, 0,
0, 0, 0, 0, 0), X3 = c(0, 0, 0, 0, 0, 0, 0, 0), X4 = c(0, 0,
0, 0, 0, 0, 0, 0), X5 = c(0, 0, 0, 0, 0, 0, 0, 0), X6 = c(0,
0, 0, 0, 0, 0, 0, 0), X7 = c(0, 0, 0, 0, 0, 0, 0, 0), X8 = c(0,
0, 0, 0, 0, 0, 1, 0), X9 = c(1, 0, 0, 0, 0, 0, 1, 1), X10 = c(1,
1, 1, 1, 1, 1, 0, 1), X11 = c(1, 1, 1, 1, 1, 1, 1, 1), X12 = c(1,
1, 1, 1, 1, 1, 1, 1)), row.names = c(NA, -8L), class = c("tbl_df",
"tbl", "data.frame"))
The first step is to calculate an infection rate for each column using a cumulative sum of positives and negatives, with the lowest value for the transition from positive to negative. To do this I am using the function rcumsum from the packageFSA
The equation is:
Infection rate = cumulative positive/(cumulative positive + cumulative negative)
and I have (a bit verbose, it definitely can be shortened I think):
Scored_Data <- as.tibble(Scored_Data)
Neg_Col <- colSums(Scored_Data)
Pos_Col <- colSums(Scored_Data == 0)
Cum_Neg_Col <- cumsum(Neg_Col)
Cum_Pos_Col <- rcumsum(Pos_Col)
Inf_Rate <- Cum_Pos_Col/(Cum_Pos_Col + Cum_Neg_Col)
Which gives me
Dput(Inf_Rate)
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 1.00000000 0.92857143 0.60000000 0.08333333 0.00000000 0.00000000
Now the difficulty comes in calculating the next step. For this I need to find the values that are immediately above and below 50% to calculate the proportional distance (PD). The equation is:
PD = (Value above 0.5 - 0.5)/(Value above 0.5 - Value below 0.5)
In my data the value above 50% is X9 and the value below is X10
I have:
PD <- for (i in Inf_Rate){
for (j in Inf_Rate){
if_else(i >= 0.5 & j < 0.5, (i-0.5)/(i-j), NULL)
}
}
Which (to my naive, non-mathematical mind) should work like this... look through the data and find the values that meet these conditions, return the PD or nothing. But this only returns "NULL".
I'm sure there is a way to get this to work, but what is it?