<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/determine-the-cv-of-a-calculated-lab-reportable-bioavailable-testosterone/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://labrtorian.com/2017/08/07/determine-the-cv-of-a-calculated-lab-reportable-bioavailable-testosterone/"> The Lab-R-torian</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Background</h2>
<p>At the <a rel="nofollow" href="http://www.aacc.org/" target="_blank">AACC</a> meeting last week, some of my friends were bugging me that I had not made a blog post in 10 months. Without getting into it too much, let's just say I can blame Cerner. Thanks also to a prod from <a rel="nofollow" href="https://www.linkedin.com/in/brian-netzel-108687105" target="_blank">a friend</a>, here is an approach to a fairly common problem.</p>
<p>We all report calculated quantities out of our laboratories–quantities such as LDL cholesterol, non-HDL cholesterol, aldosterone:renin ratio, free testosterone, eGFR etc. How does one determine the precision (i.e. imprecision) of a calculated quantity. While earlier in my life, I might go to the trouble of trying to <a rel="nofollow" href="https://doi.org/10.1016/j.clinbiochem.2006.12.014" target="_blank">do such calculations analytically</a> using the rules of <a rel="nofollow" href="https://en.wikipedia.org/wiki/Propagation_of_uncertainty" target="_blank">error propagation</a>, in my later years, I am more pragmatic and I'm happy to use a computational approach.</p>
<p>In this example, we will model the precision in calculated bioavailable testosterone (CBAT). Without explanation, I provide an R function for CBAT (and free testosterone) where testosterone is reported in nmol/L, sex hormone binding globulin (SHBG) is reported in nmol/L, and albumin is reported in g/L. Using the Vermeulen Equation as discussed in <a rel="nofollow" href="https://doi.org/10.1016/S0009-9120(03)00076-6" target="_blank">this publication</a>, you can calculate CBAT as follows:</p>

<pre>cbat &lt;- function(TT,SHBG,ALB = 43){
    Kalb &lt;- 3.6*10^4
    Kshbg &lt;- 10^9
    N &lt;- 1 + Kalb*ALB/69000
    a &lt;- N*Kshbg
    b &lt;- N + Kshbg*(SHBG - TT)/10^9
    c &lt;- -TT/10^9
    FT &lt;- (-b + sqrt(b^2 - 4*a*c))/(2*a)*10^9
    cbat &lt;- N*FT
    return(list(free.T = FT, cbat = cbat))
}</pre>

<p>To sanity-check this, we can use <a rel="nofollow" href="http://www.issam.ch/freetesto.htm" target="_blank">this online calculator</a>. Taking a typical male testosterone of 20 nmol/L, an SHBG of 50 nmol/L and an albumin of 43 g/L, we get the following:</p>

<pre>cbat(20,50)</pre>

<pre>## $free.T
## [1] 0.3273049
## 
## $cbat
## [1] 7.670319</pre>

<p>which is confirmed by the online calculator. Because the function is vectorized, we an submit a vector of testosterone results and SHBG results and get a vector of CBAT results.</p>

<pre>cbat(c(10,20,30), c(40,50,60))</pre>

<pre>## $free.T
## [1] 0.1738837 0.3273049 0.4661380
## 
## $cbat
## [1]  4.074926  7.670319 10.923842</pre>

<h2>Precision of Components</h2>
<p>We now need some precision data for the three components. However, in our lab, we just substitute 43 g/L for the albumin, so we will leave that term out of the analysis and limit our precision calculation to testosterone and SHBG. This will allow us to present the precision as surface plots as a function of total testosterone and SHBG.</p>
<p>We do testosterone by LC-MS/MS using <a rel="nofollow" href="https://doi.org/10.1016/j.cca.2012.10.007" target="_blank">Deborah French's method</a>. In the last three months, the precision has been 3.9% at 0.78 nmol/L, 5.5% at 6.7 nmol/L, 5.2% at 18.0 nmol/L, and 6.0% at 28.2 nmol/L. We are using the Roche Cobas e601 SHBG method which, according to the package insert, has precision of 1.8% at 14.9 nmol/L, 2.1 % at 45.7 nmol/L, and 4.0% at 219 nmol/L. </p>

<pre>cv.tt &lt;- c(3.9, 5.5, 5.2, 6.0)
conc.tt &lt;- c(0.78, 6.7, 18.0, 28.2)
tt.df &lt;- data.frame(conc.tt,cv.tt)

plot(cv.tt ~ conc.tt, data = tt.df,
                    main = &quot;Precision Profile of Testosterone&quot;,
                    xlab = &quot;Testosterone (nmol/L)&quot;,
                    ylab = &quot;CV Testosterone (%)&quot;,
                    ylim = c(0,8),
                    type = &quot;o&quot;)</pre>

<p><img alt="plot of chunk unnamed-chunk-4" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-4-1.png?w=604" data-recalc-dims="1"></p>

<pre>cv.shbg &lt;- c(1.8, 2.1, 4.0)
conc.shbg &lt;- c(14.9,45.7,219)
shbg.df &lt;- data.frame(cv.shbg, conc.shbg)
plot(cv.shbg ~ conc.shbg, data = shbg.df,
                    main = &quot;Precision Profile of SHBG&quot;,
                    xlab = &quot;SHBG (nmol/L)&quot;,
                    ylab = &quot;CV SHGB (%)&quot;,
                    ylim = c(0,5),
                    type = &quot;o&quot;)</pre>

<p><img alt="plot of chunk unnamed-chunk-4" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-4-2.png?w=604" data-recalc-dims="1"></p>
<h2>Build Approximation Functions</h2>
<p>We will want to generate linear interpolations of these precision profiles. Generally, we might watnt to use non-linear regression to do this but I will just linearly interpolate with the <code>approxfun()</code> function. This will allow us to just call a function to get the approximate CV at concentrations other than those for which we have data.</p>

<pre>tt.fun &lt;- approxfun(x = tt.df$conc.tt, y = tt.df$cv.tt)
shbg.fun &lt;- approxfun(x = shbg.df$conc.shbg, y = shbg.df$cv.shbg)</pre>

<p>Now, if we want to know the precision of SHBG at, say, 100 nmol/L, we can just write,</p>

<pre>shbg.fun(100)</pre>

<pre>## [1] 2.695326</pre>

<p>to obtain our precision result.</p>
<h2>Random Simulation</h2>
<p>Now let's build a grid of SHBG and total testosterone (TT) values at which we will calculate the precision for CBAT. </p>

<pre>shbg &lt;- seq(from = 15, to = 200, by = 5)
tt &lt;- seq(from = 1, to = 28, by = 1)</pre>

<p>At each point on the grid, we will have to generate, say, 100000 random TT values and 100000 random SHBG values with the appropriate precision and then calculate the expected precision of CBAT at those concentrations.</p>
<p>Let's do this for a single pair of concentrations by way of example modelling the random analytical error as Gaussian using the <code>rnorm()</code> function.</p>

<pre># [SHBG] = 15 nmol/L
# [TT] = 5.0 nmol/L
set.seed(100) #just to get consistent results
rng.tt &lt;- rnorm(100000, mean = 5.0, sd = tt.fun(5.0)/100*5.0)
rng.shbg &lt;- rnorm(100000, mean = 15, sd = shbg.fun(15)/100*15)
rng.cbat &lt;- cbat(rng.tt, rng.shbg)
cv.cbat &lt;- sd(rng.cbat$cbat)/mean(rng.cbat$cbat)*100
cv.cbat</pre>

<pre>## [1] 5.30598</pre>

<p>So, we can build the process of calculating the CV of CBAT into a function as follows:</p>

<pre>cbat.cv &lt;- function(TT, SHBG, N = 100000){
  rng.tt &lt;- rnorm(N, mean = TT, sd = tt.fun(TT)/100*TT)
  rng.shbg &lt;- rnorm(N, mean = SHBG, sd = shbg.fun(SHBG)/100*SHBG)
  rng.cbat &lt;- cbat(rng.tt, rng.shbg)
  cv &lt;- sd(rng.cbat$cbat)/mean(rng.cbat$cbat)*100
  return(cv)
}</pre>

<p>Now, we can make a matrix of the data for presenting a plot, calculating the CV and appending it to the dataframe.</p>

<pre>cv.grid &lt;- expand.grid(tt, shbg)
names(cv.grid) &lt;- c(&quot;tt&quot;, &quot;shbg&quot;)
cv.grid$cv.cbat &lt;- mapply(cbat.cv, cv.grid$tt, cv.grid$shbg)</pre>

<p>Now make plot using the <code>wireframe()</code> function.</p>

<pre>library(lattice)
wireframe(cv.cbat ~ tt*shbg, data = cv.grid,
          xlab = &quot;Testo \n (nmol/L)&quot;,
          ylab = &quot;SHBG \n (nmol/L)&quot;,
          zlab = &quot;CV \n (%)&quot;,
          drape = TRUE,
          colorkey = TRUE,
          col.regions = colorRampPalette(c(&quot;blue&quot;, &quot;red&quot;, &quot;yellow&quot;))(100),
          scales = list(arrows=FALSE,cex=.5,tick.number = 10)
          )</pre>

<p><img alt="plot of chunk unnamed-chunk-11" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-11-1.png?w=604" data-recalc-dims="1"></p>
<p>This shows us that the CV of CBAT ranges from about 4–8% over the TT and SHBG ranges we have looked at. </p>
<h2>Conclusion</h2>
<p>We have determined the CV of calculated bioavailable testosterone using random number simulations using empirical CV data and produced a surface plot of CV. This allows us to comment on the CV of this lab reportable as a function of the two variables by which it is determined.</p>

<p><strong>Parting Thought on Monte Carlo Simulations</strong></p>
<p><em>The die is cast into the lap, but its every decision is from the LORD.</em></p>
<p>(Prov 16:33)</p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://labrtorian.com/2017/08/07/determine-the-cv-of-a-calculated-lab-reportable-bioavailable-testosterone/"> The Lab-R-torian</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/qHe6UAsvGxQ" height="1"></a></div>