Predicting when an output might happen in time in R

Question

I have a dataset of peoples blood results in R:

        ID Result date    
        A1 80     01/01/2006
        A1 70     01/01/2009
        A1 61     01/01/2010
        A1 30     01/01/2018
        A1 28     01/01/2022
        B2 40     01/01/2006
        B2 30     01/01/2012
        B2 25     01/01/2015
        B2 10     01/01/2020
        G3 28     01/01/2009
        G3 27     01/01/2014
        G3 25     01/01/2013
        G3 24     01/01/2011
        G3 22     01/01/2019

I would like to plot these and then do a linear regression and find the predicted date at which their blood result would hit 15 for each person. In reality there are multiple blood tests per year so I would try and get the median result per year and then plot those and get the predicted year they hit 15.

I can do the median per year and use the lm function along the lines of:

filtered$date_of_bloods <-format(filtered$date_of_bloods,format="%Y")
#split into individual ID groups
a <- with(filtered, split(filtered, list(ID)))

#aggregate median results per year 
medianfunc <- function(y) {aggregate(results ~ date_of_bloods, data = y, median)}
medians <- sapply(a, medianfunc)

# do lm per ID cohort and get slope of lines 
g<- as.data.frame(medians)
coefLM <- function(x) {coef(lm(date_of_bloods ~ results, data = x))}
coefs<- sapply(g, coefLM)

But I am unsure who to extract and get the prediction for when the slope of the line would hit 15.

Any help would be much appreciated

score 1 · Answer 1 · answered Sep 12 '22 at 16:24

Here is solution relying on dplyr, splitting the data into lists to feed to lapply, where you can define a custom function. For this simple linear regression, we use the coefficients to calculate X for when Y = 15.

You should round the years to your desired level of precision. Note that this does NOT check the assumptions for linear regression, or report on the goodness-of-fit for your given models.

# Use example data from OP
df = structure(list(ID = c("A1", "A1", "A1", "A1", "A1", "B2", 
"B2", "B2", "B2", "G3", "G3", "G3", "G3", "G3"), Result = c(80L, 70L, 
61L, 30L, 28L, 40L, 30L, 25L, 10L, 28L, 27L, 25L, 24L, 22L), 
date = c("01/01/2006", "01/01/2009", "01/01/2010", "01/01/2018", 
"01/01/2022", "01/01/2006", "01/01/2012", "01/01/2015", "01/01/2020", 
"01/01/2009", "01/01/2014", "01/01/2013", "01/01/2011", "01/01/2019"
)), class = "data.frame", row.names = c(NA, -14L))

library(dplyr)
library(lubridate)

df %>% 
  # Create a column holding year for each ID
  mutate(date = dmy(date)) %>% 
  mutate(year = year(date)) %>% 
  # Make yearly concentrations for each id by year, using e.g. median
  group_by(ID, year) %>% 
  summarise(conc = median(Result), .keep_all = TRUE) %>% 
  # Split each individual into separate lists
  split(f = .$ID) %>% 
  lapply(., function(x) {
    
    # Create LM object
    mod = lm(conc ~ year, data = x)
    
    # Use coefficients to find when y = 15
    when_is_15 = (15 - coef(mod)[1]) / coef(mod)[2]
    
    # Wrap up the results
    data.frame(ID = unique(x$ID),
               pred_year = as.numeric(when_is_15))
    
  }) %>% 
  # Combine all the lists into a single data frame
  bind_rows()

#>   ID pred_year
#> 1 A1  2024.246
#> 2 B2  2018.595
#> 3 G3  2035.313

^{Created on 2022-09-12 with reprex v2.0.2}

@tacrolimus, this is handled by taking the median of those samples. You can change that by editing the line `summarise(conc = median(Result), .keep_all = TRUE) `, to whichever summary statistic best fits your work — mhovd, Sep 13 '22 at 10:47
Many thanks for taking the time to answer. This worked really well. I have selected the other answer as it uses base R and I work within an airlocked environment making it difficult to update and import packages. But this was great thank you! — tacrolimus, Sep 13 '22 at 10:47

score 1 · Accepted Answer · answered Sep 12 '22 at 16:36

Here's a base R answer:

filtered <- read.table(text = "ID results date_of_bloods    
                A1 80     01/01/2006
                A1 70     01/01/2009
                A1 61     01/01/2010
                A1 30     01/01/2018
                A1 28     01/01/2022
                B2 40     01/01/2006
                B2 30     01/01/2012
                B2 25     01/01/2015
                B2 10     01/01/2020
                G3 28     01/01/2009
                G3 27     01/01/2014
                G3 25     01/01/2013
                G3 24     01/01/2011
                G3 22     01/01/2019", header = TRUE)
filtered$date_of_bloods <- as.Date(filtered$date_of_bloods, '%m/%d/%Y')

filtered$date_of_bloods <- format(filtered$date_of_bloods,format="%Y") |> 
  as.integer()
#split into individual ID groups
a <- with(filtered, split(filtered, list(ID)))

#aggregate median results per year 
medianfunc <- function(y) {aggregate(results ~ date_of_bloods, data = y, median)}
medians <- sapply(a, medianfunc)

# do lm per ID cohort and get slope of lines 
g <- as.data.frame(medians)
coefLM <- function(x) {coef(lm(results ~ date_of_bloods, data = x))}
coefs<- sapply(g, coefLM)  

coefs <- 
  rbind(coefs, 
        year_for_15 = (15 - coefs["(Intercept)", ]) / coefs["date_of_bloods", ])

Giving:

coefs
                     A1          B2           G3
(Intercept)    6998.650 4263.381995  953.8239437
date_of_bloods   -3.450   -2.104623   -0.4612676
year_for_15    2024.246 2018.595376 2035.3129771

I changed your code for the year to make it an integer, that way the coefLM function can work in the normal dependent variable ~ independent variable order. Then it's just y = mx + b implies x0 = (15 - b) / m.

Many thanks for taking the time to answer. This worked really well and makes clear sense to me! — tacrolimus, Sep 13 '22 at 10:39

Predicting when an output might happen in time in R

2 Answers2

Linked