How to create a (total) Ghrelin calibration curve?
- Version 3
- by (unknown)
- Version 23
- by (unknown)
Deletions or items before changed
Additions or items after changed
1 | - | + | == Computing Ghrelin calibration curve using R ==
|
|
---|---|---|---|---|
2 | ||||
3 | + | ===This article describes creating calibration curve, and measuring the total-ghrelin concentration ===
|
||
4 | ||||
5 | + | 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.
|
||
6 | ||||
7 | - | < |
+ | We will assume a 4-parameter log-logistic model:
|
8 | - | + | ||
9 | - | + | [[Image(LL4.png, 300px)]]
|
|
10 | + | |||
11 | + | The R code is attached below.
|
||
12 | + | |||
13 | + | The example assumes the data to be available in file "ghrelin_conc std_a std_b avg.csv"
|
||
14 | + | |||
15 | + | The measured data:
|
||
16 | + | |||
17 | + | || '''Ghrelin (ng/ml)''' || '''Standard a''' || '''Standard b''' ||
|
||
18 | + | || 1000000 ||-0.040596823||-0.052699697 ||
|
||
19 | + | || 100000 ||0.136105144||0.119766263 ||
|
||
20 | + | || 10000 ||0.61356354||0.606906959 ||
|
||
21 | + | || 1000 ||0.846543873||0.839887292 ||
|
||
22 | + | || 100 ||0.887693646||0.88345764 ||
|
||
23 | + | ||0||0.896770802||0.896165658 ||
|
||
24 | + | |||
25 | + | {{{
|
||
26 | + | |||
27 | + | ##### Install libraries
|
||
28 | + | install.packages("drc")
|
||
29 | + | install.packages("sfsmisc")
|
||
30 | + | require(drc)
|
||
31 | + | library(sfsmisc)
|
||
32 | + | |||
33 | + | ##### Read the data
|
||
34 | + | hormone.data <- read.csv("ghrelin_conc std_a std_b avg.csv")
|
||
35 | + | hormone.data <- hormone.data[,1:3]
|
||
36 | + | colnames(hormone.data)[1:3] <- c("Concentration","Response_1", "Response_2")
|
||
37 | + | |||
38 | + | ##### Reorganize the data
|
||
39 | + | hormone.data <- reshape(hormone.data, varying=c("Response_1","Response_2"), direction="long", v.names=c("Response"))
|
||
40 | + | hormone.data <- hormone.data[,c("Concentration", "Response")]
|
||
41 | + | |||
42 | + | ##### Fitting the model (4-parameter log-logistic function)
|
||
43 | + | hormone.data.model <- drm(Response ~ Concentration, data = hormone.data, fct = LL.4())
|
||
44 | + | summary(hormone.data.model)
|
||
45 | + | |||
46 | + | }}}
|
||
47 | + | |||
48 | + | The resultant parameters of a log-logistic equation are:
|
||
49 | + | |||
50 | + | {{{
|
||
51 | + | Model fitted: Log-logistic (ED50 as parameter) (4 parms)
|
||
52 | + | Parameter estimates:
|
||
53 | + | Estimate Std. Error t-value p-value
|
||
54 | + | b:(Intercept) 9.5057e-01 2.2294e-02 4.2638e+01 0
|
||
55 | + | c:(Intercept) -7.6010e-02 6.9075e-03 -1.1004e+01 0
|
||
56 | + | d:(Intercept) 8.9163e-01 3.3216e-03 2.6843e+02 0
|
||
57 | + | e:(Intercept) 2.5221e+04 7.7727e+02 3.2448e+01 0
|
||
58 | + | }}}
|
||
59 | + | |||
60 | + | The calibration curve can be plotted using the commands below:
|
||
61 | + | |||
62 | + | {{{
|
||
63 | + | ##### Plotting a nice plot
|
||
64 | + | par(pty="s", mar=c(5,5,1,1))
|
||
65 | + | plot(hormone.data.model, type="confidence", cex.lab=2, axes=F, xlim=c(-10,10^6))
|
||
66 | + | axis(side=1, at=hormone.data[1:6,1], labels=pretty10exp(hormone.data[1:6,1]), cex.axis=1.2)
|
||
67 | + | axis(side=2, at=seq(0,1,0.2), labels=seq(0,1,0.2))
|
||
68 | + | plot(hormone.data.model, type="all", add=T, pch=21, col="red", lwd=1, cex=2, bg="green")
|
||
69 | + | }}}
|
||
70 | + | |||
71 | + | [[Image(Ghrelin.png)]]
|
||
72 | + | |||
73 | + | 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.
|
||
74 | + | |||
75 | + | [[Image(LL4-inv.png, 230px)]]
|
||
76 | + | |||
77 | + | 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.
|
||
78 | + | |||
79 | + | {{{
|
||
80 | + | ##### Computing the concentration from the response, for instance for a response=0.1, and alpha=1-0.95
|
||
81 | + | ED(hormone.data.model, respLev=0.1, interval="delta", type="absolute", level=0.95)
|
||
82 | + | }}} |