---
title: "Non-linear regression of Elisa data"
author: "Burnet Bioinformatics Group"
date: "`r Sys.Date()`"
output:
  html_document:
    theme: simplex
    toc: true
    code_folding: hide
    fig_width: 7
    fig_height: 7
    toc_float:
      collapsed: true
---

## Introduction

Need more accurate approaches to quantify differences in plasma concentration of
analytes that are tolerant of outliers.
My suggestion was to generate non-linear regression curves of Elisa data from each
visit and use a statistical test to determine the difference.


```{r,libs}

library("kableExtra")

```

## Example regression

Exponential decay non-linear regression in R

```{r,example1}

set.seed(1)
x <- 1:20
y <- 5 * exp(-0.3 * x) + 0.5 + rnorm(20, 0, 0.2)
df <- data.frame(x, y)

df |>
  kbl(caption="Synthetic dataset") |>
  kable_paper("hover", full_width = F)

fit <- nls(
  y ~ a * exp(-b * x) + c,
  data = df,
  start = list(a = 5, b = 0.3, c = 0.1)
)
summary(fit)

plot(df$x, df$y)
curve(coef(fit)["a"] * exp(-coef(fit)["b"] * x) + coef(fit)["c"],
      add = TRUE, lwd = 2, col = "red")

```

That looks great.
Now try fitting a combined model and use an ANOVA (nested-model comparison).

```{r,combinedmodel_example}

set.seed(100)
x <- 1:20
y <-c( 5 * exp(-0.3 * x) + 0.5 + rnorm(20, 0, 0.2) ,
     10 * exp(-0.5 * x) + 0.5 + rnorm(20, 1, 0.2) )
x <- rep(x,2)

df0 <- data.frame(x, y)
df0$group <- c(rep("A",20),rep("B",20))

df0 |>
  kbl(caption="Synthetic dataset2") |>
  kable_paper("hover", full_width = F)

fit_null <- nls(
  y ~ a * exp(-b * x) + c,
  data = df0,
  start = list(a=5, b=0.3, c=0.5)
)

fit_alt <- nls(
  y ~ (a1 + (a2 - a1) * (group == "B")) *
       exp(-(b1 + (b2 - b1) * (group == "B")) * x) +
       (c1 + (c2 - c1) * (group == "B")),
  data = df0,
  start = list(a1=5, a2=5, b1=0.3, b2=0.3, c1=0.5, c2=0.5)
)

anova(fit_null, fit_alt)

```

## Data for patient 82

Now we will look at some real data.

```{r,read1}

df1 <- read.table("pat82_visit1.tsv",header=TRUE)
x1 <- log10(rep(df1$plasma_dilution,3))
y1 <- c(df1$rep1,df1$rep2,df1$rep3)
df1 <- data.frame(x1, y1)
df1$group <- "A"

df1 |>
  kbl(caption="Patient 82 visit 1") |>
  kable_paper("hover", full_width = F)

fit1 <- nls(
  y1 ~ a * exp(-b * x1) + c,
  data = df1,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit1)

df2 <- read.table("pat82_visit2.tsv",header=TRUE)
x2 <- log10(rep(df2$plasma_dilution,3))
y2 <- c(df2$rep1,df2$rep2,df2$rep3)
df2 <- data.frame(x2, y2)
df2$group <- "B"

df2 |>
  kbl(caption="Patient 82 visit 2") |>
  kable_paper("hover", full_width = F)

fit2 <- nls(
  y2 ~ a * exp(-b * x2) + c,
  data = df1,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit2)

plot(df1$x,df1$y,pch=15,main="all data",
  xlab="Plasma dilution factor",
  ylab="Absorbance")

points(df2$x,df2$y,col="blue",pch=17)

curve(coef(fit1)["a"] * exp(-coef(fit1)["b"] * x) + coef(fit1)["c"],
      add = TRUE, lwd = 2, col = "black", lty=2)

curve(coef(fit2)["a"] * exp(-coef(fit2)["b"] * x) + coef(fit2)["c"],
      add = TRUE, lwd = 2, col = "blue", lty=2)

```

It looks like there are some outliers there.
Let's quantify their variation from the regression curve and exclude the bad ones.

```{r,outliers}

df1$expect1 <- coef(fit1)["a"] * exp(-coef(fit1)["b"] * df1$x1) + coef(fit1)["c"]
df1$delta1 <- abs(df1$y1 - df1$expect1)
df1$cv1 <- abs(df1$delta / df1$expect1)
df1clean <- subset(df1,cv1<3)

df1clean |>
  kbl(caption="Patient 82 visit 1 clean") |>
  kable_paper("hover", full_width = F)

df2$expect2 <- coef(fit2)["a"] * exp(-coef(fit2)["b"] * df2$x2) + coef(fit2)["c"]
df2$delta1 <- abs(df2$y2 - df2$expect2)
df2$cv2 <- abs(df2$delta / df2$expect2)
df2clean <- subset(df2,cv2<3)

df2clean |>
  kbl(caption="Patient 82 visit 2 clean") |>
  kable_paper("hover", full_width = F)

```

Now that we have excluded the outliers, we can run regression again.

```{r,fit2}

fit1c <- nls(
  y1 ~ a * exp(-b * x1) + c,
  data = df1clean,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit1c)

fit2c <- nls(
  y2 ~ a * exp(-b * x2) + c,
  data = df2clean,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit2c)

plot(df1clean$x,df1clean$y,pch=15, main="outliers removed",
  xlab="Plasma dilution factor",
  ylab="Absorbance")

points(df2clean$x,df2clean$y,col="blue",pch=17)

curve(coef(fit1c)["a"] * exp(-coef(fit1c)["b"] * x) + coef(fit1c)["c"],
      add = TRUE, lwd = 2, col = "black", lty=2)

curve(coef(fit2c)["a"] * exp(-coef(fit2c)["b"] * x) + coef(fit2c)["c"],
      add = TRUE, lwd = 2, col = "blue", lty=2)

```

This is looking a lot better.


Calculate difference between curves.

Predict difference between curves.
Reliable range is 2-2.5.
Calculate foldchange.

```{r,predict82}

newx <- data.frame(x1 = seq(2, 2.5, length.out = 5))
predicted1 <- predict(fit1c, newdata = newx)
newx <- data.frame(x2 = seq(2, 2.5, length.out = 5))
predicted2 <- predict(fit2c, newdata = newx)
pred_diff <- (predicted2-predicted1)/predicted1
pred_diff
FOLDCHANGE=signif(mean(pred_diff),3)
message(paste("foldchange is",FOLDCHANGE))

```

## Patient 82 combined model with ANOVA model comparison

```{r,pat82_combined}

df1clean$group <- "A"
df2clean$group <- "B"

colnames(df2clean) <- colnames(df1clean)
df0 <- rbind(df1clean,df2clean)

df0 |>
  kbl(caption="Patient 82 visit 1 and 2 clean") |>
  kable_paper("hover", full_width = F)

fit_null <- nls(
  y1 ~ a * exp(-b * x1) + c,
  data = df0,
  start = list(a=1e6, b=0.3, c=5e4)
)
fit_null

fit_alt <- nls(
  y1 ~ (a1 + (a2 - a1) * (group == "B")) *
       exp(-(b1 + (b2 - b1) * (group == "B")) * x1) +
       (c1 + (c2 - c1) * (group == "B")),
  data = df0,
  start = list(a1=1e6, a2=1e6, b1=0.3, b2=0.3, c1=5e4, c2=5e4)
)
fit_alt

anova(fit_null, fit_alt)

```

Differences in these curves is not significant (p=0.95).

## Data for patient 440

```{r,read2}

df1 <- read.table("pat440_visit1.tsv",header=TRUE)
x1 <- log10(rep(df1$plasma_dilution,3))
y1 <- c(df1$rep1,df1$rep2,df1$rep3)
df1 <- data.frame(x1, y1)

df1 |>
  kbl(caption="Patient 440 visit 1") |>
  kable_paper("hover", full_width = F)

fit1 <- nls(
  y1 ~ a * exp(-b * x1) + c,
  data = df1,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit1)

df2 <- read.table("pat440_visit2.tsv",header=TRUE)
x2 <- log10(rep(df2$plasma_dilution,3))
y2 <- c(df2$rep1,df2$rep2,df2$rep3)
df2 <- data.frame(x2, y2)

df2 |>
  kbl(caption="Patient 440 visit 2") |>
  kable_paper("hover", full_width = F)

fit2 <- nls(
  y2 ~ a * exp(-b * x2) + c,
  data = df1,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit2)

MAX=max(c(df1$y1,df2$y2))

plot(df1$x,df1$y,pch=15,main="all data", ylim=c(0,MAX),
  xlab="Plasma dilution factor",
  ylab="Absorbance")

points(df2$x,df2$y,col="blue",pch=17)

curve(coef(fit1)["a"] * exp(-coef(fit1)["b"] * x) + coef(fit1)["c"],
      add = TRUE, lwd = 2, col = "black", lty=2)

curve(coef(fit2)["a"] * exp(-coef(fit2)["b"] * x) + coef(fit2)["c"],
      add = TRUE, lwd = 2, col = "blue", lty=2)

```

It looks like there are some outliers there.
Let's quantify their variation from the regression curve and exclude the bad ones.

```{r,outliers2}
df1$expect1 <- coef(fit1)["a"] * exp(-coef(fit1)["b"] * df1$x1) + coef(fit1)["c"]
df1$delta1 <- abs(df1$y1 - df1$expect1)
df1$cv1 <- abs(df1$delta / df1$expect1)
df1clean <- subset(df1,cv1<3)

df1clean |>
  kbl(caption="Patient 440 visit 1 clean") |>
  kable_paper("hover", full_width = F)

df2$expect2 <- coef(fit2)["a"] * exp(-coef(fit2)["b"] * df2$x2) + coef(fit2)["c"]
df2$delta1 <- abs(df2$y2 - df2$expect2)
df2$cv2 <- abs(df2$delta / df2$expect2)
df2clean <- subset(df2,cv2<3)

df2clean |>
  kbl(caption="Patient 440 visit 2 clean") |>
  kable_paper("hover", full_width = F)

```

Now that we have excluded the outliers, we can run regression again.

```{r,fit3}

fit1c <- nls(
  y1 ~ a * exp(-b * x1) + c,
  data = df1clean,
  start = list(a = 1e6, b = 0.3, c = 1e5)
)
summary(fit1c)

fit2c <- nls(
  y2 ~ a * exp(-b * x2) + c,
  data = df2clean,
  start = list(a = 5e6, b = 0.3, c = 1e5)
)
summary(fit2c)

MAX=max(c(df1clean$y1,df2clean$y2))

plot(df1clean$x,df1clean$y,pch=15, main="outliers removed",ylim=c(0,MAX),
  xlab="Plasma dilution factor",
  ylab="Absorbance")

points(df2clean$x,df2clean$y,col="blue",pch=17)

curve(coef(fit1c)["a"] * exp(-coef(fit1c)["b"] * x) + coef(fit1c)["c"],
      add = TRUE, lwd = 2, col = "black", lty=2)

curve(coef(fit2c)["a"] * exp(-coef(fit2c)["b"] * x) + coef(fit2c)["c"],
      add = TRUE, lwd = 2, col = "blue", lty=2)

```

These curves are different, but we need to model and test them.

This is looking a lot better.

## Patient 440 combined model with ANOVA model comparison

```{r,pat440_combined}

df1clean$group <- "A"
df2clean$group <- "B"

colnames(df2clean) <- colnames(df1clean)
df0 <- rbind(df1clean,df2clean)

df0 |>
  kbl(caption="Patient 440 visit 1 and 2") |>
  kable_paper("hover", full_width = F)

fit_null <- nls(
  y1 ~ a * exp(-b * x1) + c,
  data = df0,
  start = list(a=3e6, b=0.3, c=1e5)
)
fit_null

fit_alt <- nls(
  y1 ~ (a1 + (a2 - a1) * (group == "B")) *
       exp(-(b1 + (b2 - b1) * (group == "B")) * x1) +
       (c1 + (c2 - c1) * (group == "B")),
  data = df0,
  start = list(a1=3e6, a2=3e6, b1=0.3, b2=0.3, c1=1e5, c2=1e5)
)
fit_alt

anova(fit_null, fit_alt)

anova(fit_null, fit_alt)$Pr[2]
```

The difference between curves is significant (p=9.2e-22).

Predict difference between curves.
Reliable range is 2-2.5.
Calculate foldchange.

```{r,predict440}

newx <- data.frame(x1 = seq(2, 2.5, length.out = 5))
predicted1 <- predict(fit1c, newdata = newx)
newx <- data.frame(x2 = seq(2, 2.5, length.out = 5))
predicted2 <- predict(fit2c, newdata = newx)
pred_diff <- (predicted2-predicted1)/predicted1
pred_diff
FOLDCHANGE=signif(mean(pred_diff),3)
message(paste("foldchange is",FOLDCHANGE))

```
