Computing Ghrelin calibration curve using R
The calibration curve can be created using R, and libracy “drc”. Optionally, one can use library “sfmisc” for formatting of the labels on plot axis.
We will assume a 4-parameter log-logistic model:
The R code is attached below.
The example assumes the data to be available in file “ghrelin_conc std_a std_b avg.csv”
The measured data:
| Ghrelin (ng/ml) | Standard a | Standard b |
| 1000000 | -0.040596823 | -0.052699697 |
| 100000 | 0.136105144 | 0.119766263 |
| 10000 | 0.61356354 | 0.606906959 |
| 1000 | 0.846543873 | 0.839887292 |
| 100 | 0.887693646 | 0.88345764 |
| 0 | 0.896770802 | 0.896165658 |
##### Install libraries
install.packages("drc")
install.packages("sfsmisc")
require(drc)
library(sfsmisc)
##### Read the data
hormone.data <- read.csv("ghrelin_conc std_a std_b avg.csv")
hormone.data <- hormone.data[,1:3]
colnames(hormone.data)[1:3] <- c("Concentration","Response_1", "Response_2")
##### Reorganize the data
hormone.data <- reshape(hormone.data, varying=c("Response_1","Response_2"), direction="long", v.names=c("Response"))
hormone.data <- hormone.data[,c("Concentration", "Response")]
##### Fitting the model (4-parameter log-logistic function)
hormone.data.model <- drm(Response ~ Concentration, data = hormone.data, fct = LL.4())
summary(hormone.data.model)
The resultant parameters of a log-logistic equation are:
Model fitted: Log-logistic (ED50 as parameter) (4 parms)
Parameter estimates:
Estimate Std. Error t-value p-value
b:(Intercept) 9.5057e-01 2.2294e-02 4.2638e+01 0
c:(Intercept) -7.6010e-02 6.9075e-03 -1.1004e+01 0
d:(Intercept) 8.9163e-01 3.3216e-03 2.6843e+02 0
e:(Intercept) 2.5221e+04 7.7727e+02 3.2448e+01 0
The calibration curve can be plotted using the commands below:
##### Plotting a nice plot par(pty="s", mar=c(5,5,1,1)) plot(hormone.data.model, type="confidence", cex.lab=2, axes=F, xlim=c(-10,10^6)) axis(side=1, at=hormone.data[1:6,1], labels=pretty10exp(hormone.data[1:6,1]), cex.axis=1.2) axis(side=2, at=seq(0,1,0.2), labels=seq(0,1,0.2)) plot(hormone.data.model, type="all", add=T, pch=21, col="red", lwd=1, cex=2, bg="green")
The parameters of the eqution can be plugged into the formula below (an inverse of the model), and used in Excel, or other spreadsheet program.
However, the concentration can be also easily estimated in R using “ED” function of the “drc” library. The code below demonstrates the concentration estimated from the response of 0.1, assuming alpha=0.05. The code returns the estimation, the error, and the condfidence interval.
##### Computing the concentration from the response, for instance for a response=0.1, and alpha=1-0.95 ED(hormone.data.model, respLev=0.1, interval="delta", type="absolute", level=0.95)