<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/non-linear-regression-application-to-monoclonal-peak-integration-in-serum-protein-electrophoresis/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://labrtorian.com/2017/08/28/non-linear-regression-application-to-monoclonal-peak-integration-in-serum-protein-electrophoresis/"> 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 recently, there was an enthusiastic discussion of standardization of reporting for serum protein electrophoresis (SPEP) presented by a working group headed up by <a rel="nofollow" href="https://www.linkedin.com/in/christopher-mccudden-5136358/?ppe=1" target="_blank">Dr. Chris McCudden</a> and <a rel="nofollow" href="http://www.med.uottawa.ca/patho/eng/booth.html" target="_blank">Dr. Ron Booth</a>, both of the University of Ottawa. One of the discussions pertained to how monoclonal bands, especially small ones, should be integrated. While many use the default manual vertical gating or “drop” method offered by <a rel="nofollow" href="http://labrtorian.com/2017/08/28/non-linear-regression-application-to-monoclonal-peak-integration-in-serum-protein-electrophoresis/www.sebia.com" target="_blank">Sebia's</a> Phoresis software, <a rel="nofollow" href="https://www.pathology.med.umich.edu/faculty/dkeren" target="_blank">Dr. David Keren</a> was discussing the value of tangent skimming as a more repeatable and effective means of monoclonal protein quantitation. He was also discussing some biochemical approaches distinguishing monoclonal proteins from the background gamma proteins. </p>
<p>The drop method is essentially an eye-ball approach to where the peak starts and ends and is represented by the vertical lines and the enclosed shaded area.</p>
<p><img alt="plot of chunk unnamed-chunk-1" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-1-1.png?w=604" data-recalc-dims="1"></p>
<p>The tangent skimming approach <a rel="nofollow" href="https://www.degruyter.com/view/j/cclm.2016.54.issue-6/cclm-2015-0862/cclm-2015-0862.xml" target="_blank">is easier to make reproducible</a>. In the mass spectrometry world it is a well-developed approach with a long history and <a rel="nofollow" href="http://www.chromatographyonline.com/integration-errors-chromatographic-analysis-part-i-peaks-approximately-equal-size" target="_blank">multiple algorithms in use</a>. This is apparently <a rel="nofollow" href="https://books.google.ca/books?id=wGwoDwAAQBAJ&printsec=frontcover&dq=Chromatographic+Integration+Methods&hl=en&sa=X&ved=0ahUKEwjLpqPq1efVAhWmhlQKHcffAI8Q6AEIKDAA#v=onepage&q=Chromatographic%20Integration%20Methods&f=false" target="_blank"><em>the</em> book</a>. However, when tangent skimming is employed in SPEP, unless I am mistaken, it seems to be done by eye. The integration would look like this:</p>
<p><img alt="plot of chunk unnamed-chunk-2" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-2-1.png?w=604" data-recalc-dims="1"></p>
<p>During the discussion it was point out that peak deconvolution of the monoclonal protein from the background gamma might be preferable to either of the two described procedures. By this I mean integration as follows:</p>
<p><img alt="plot of chunk unnamed-chunk-3" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-3-1.png?w=604" data-recalc-dims="1"></p>
<p>There was discussion this procedure is challenging for number of reasons. Further, it should be noted that there will only likely be any clinical value in a deconvolution approach when the concentration of the monoclonal protein is low enough that manual integration will show poor repeatability, say &lt; 5 g/L = 0.5 g/dL. </p>
<h2>Easy Peaks</h2>
<p>Fitting samples with larger monoclonal peaks is fairly easy. Fitting tends to converge nicely and produce something meaningful. For example, using the approach I am about to show below, an electropherogram like this:</p>
<p><img alt="plot of chunk unnamed-chunk-4" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-4-1-1.png?w=604" data-recalc-dims="1"></p>
<p>with a gamma region looking like this:</p>
<p><img alt="plot of chunk unnamed-chunk-5" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-5-1.png?w=604" data-recalc-dims="1"></p>
<p>can be deconvoluted with straightforward non-linear regression (and no baseline subtraction) to yield this:</p>
<p><img alt="plot of chunk unnamed-chunk-6" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-6-1-1.png?w=604" data-recalc-dims="1"></p>
<p>and the area of the green monoclonal peak is found to be 5.3%.</p>
<h2>More Difficult Peaks</h2>
<p>What is more challenging is the problem of small monoclonals buried in normal \(\gamma\)-globulins. These could be difficult to integrate using a tangent skimming approach, particularly without image magnification. For the remainder of this post we will use a gel with a small monoclonal in the fast gamma region shown at the arrow.</p>
<p><img alt="plot of chunk unnamed-chunk-7" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-7-1-1.png?w=604" data-recalc-dims="1"></p>
<h3>Getting the Data</h3>
<p>EP data can be extracted from the PDF output from any electrophoresis software. This is not complicated and can be accomplished with <a rel="nofollow" href="http://www.cityinthesky.co.uk/opensource/pdf2svg/" target="_blank">pdf2svg</a> or <a rel="nofollow" href="https://inkscape.org/en/" target="_blank">Inkscape</a> and some Linux bash scripting. I'm sure we can get it straight from the instrument but it is not obvious to me how to do this. One could also rescan a gel and use <a rel="nofollow" href="https://imagej.nih.gov/ij/" target="_blank">ImageJ</a> to produce a densitometry scan which is discussed in the <a rel="nofollow" href="https://imagej.nih.gov/ij/docs/guide/146-30.html#infobox:Densitometry" target="_blank">ImageJ documentation</a> and on <a rel="nofollow" href="https://www.youtube.com/watch?v=JlR5v-DsTds" target="_blank">YouTube</a>. <a rel="nofollow" href="https://imagej.nih.gov/ij/developer/macro/macros.html" target="_blank">ImageJ also has a macro language</a> for situations where the same kind of processing is done repeatedly.</p>
<h3>Smoothing</h3>
<p>The data has 10284 pairs of (x,y) data. But if you blow up on it and look carefully you find that it is a series of staircases.</p>

<pre>plot(y~x, data = head(ep.data,100), type = &quot;o&quot;, cex = 0.5)</pre>

<p><img alt="plot of chunk unnamed-chunk-8" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-8-1.png?w=604" data-recalc-dims="1"></p>
<p>It turns out that this jaggedness significantly impairs attempts to numerically identify the peaks and valleys. So, I smoothed it a little using the handy <code>rle()</code> function to identify the midpoint of each step. This keeps the total area as close to its original value as possible–though this probably does not matter too much.</p>

<pre>ep.rle &lt;- rle(ep.data$y)
stair.midpoints &lt;- cumsum(ep.rle$lengths) - floor(ep.rle$lengths/2)
ep.data.sm &lt;- ep.data[stair.midpoints,]
plot(y~x, data = head(ep.data,300), type = &quot;o&quot;, cex = 0.5)
points(y~x, data = head(ep.data.sm,300), type = &quot;o&quot;, cex = 0.5, col = &quot;red&quot;)</pre>

<p><img alt="plot of chunk unnamed-chunk-9" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-9-1.png?w=604" data-recalc-dims="1"></p>
<p>Now that we are satisfied that the new data is OK, I will overwrite the original dataframe.</p>

<pre>ep.data &lt;- ep.data.sm</pre>

<h3>Transformation</h3>
<p>The units on the x and y-axes are arbitrary and come from page coordinates of the PDF. We can normalize the scan by making the x-axis go from 0 to 1 and by making the total area 1.</p>

<pre>library(Bolstad) #A package containing a function for Simpon's Rule integration
ep.data$x &lt;- ep.data$x/max(ep.data$x)
A.tot &lt;- sintegral(ep.data$x,ep.data$y)$value
ep.data$y &lt;- ep.data$y/A.tot

#sanity check
sintegral(ep.data$x,ep.data$y)$value</pre>

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

<pre>plot(y~x, data = ep.data, type = &quot;l&quot;)</pre>

<p><img alt="plot of chunk unnamed-chunk-11" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-11-1-1.png?w=604" data-recalc-dims="1"></p>
<h3>Find Extrema</h3>
<p>Using the <code>findPeaks</code> function from the <a rel="nofollow" href="https://cran.r-project.org/web/packages/quantmod/index.html" target="_blank">quantmod</a> package we can find the minima and maxima:</p>

<pre>library(quantmod)
ep.max &lt;- findPeaks(ep.data$y)
plot(y~x, data = ep.data, type = &quot;l&quot;, main = &quot;Maxima&quot;)
abline(v = ep.data$x[ep.max], col = &quot;red&quot;, lty = 2)</pre>

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

<pre>ep.min &lt;- findValleys(ep.data$y)
plot(y~x, data = ep.data, type = &quot;l&quot;, main = &quot;Minima&quot;)
abline(v = ep.data$x[ep.min], col = &quot;blue&quot;, lty = 2)</pre>

<p><img alt="plot of chunk unnamed-chunk-12" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-12-2.png?w=604" data-recalc-dims="1"></p>
<p>Not surprisingly, there are some extraneous local extrema that we do not want. I simply manually removed them. Generally, this kind of thing could be tackled with more smoothing of the data prior to analysis.</p>

<pre>ep.max &lt;- ep.max[-1]
ep.min &lt;- ep.min[-c(1,length(ep.min))]</pre>

<h2>Fitting</h2>
<p>Now it's possible with the <code>nls()</code> function to fit the entire SPEP with a series of Gaussian curves simultaneously. It works just fine (provided you have decent initial estimates of \(\mu_i\) and \(\sigma_i\)) but there is no particular clinical value to fitting the albumin, \(\alpha_1\), \(\alpha_2\), \(\beta_1\) and \(\beta_2\) domains with Gaussians. What is of interest is separately quantifying the two peaks in \(\gamma\) with two separate Gaussians so let's isolate the \(\gamma\) region based on the location of the minimum between \(\beta_2\) and \(\gamma\).</p>
<h3>Isolate the \(\gamma\) Region</h3>

<pre>gamma.ind &lt;- max(ep.min):nrow(ep.data)
gamma.data &lt;- data.frame(x = ep.data$x[gamma.ind], y = ep.data$y[gamma.ind])
plot(y ~ x, gamma.data, type  = &quot;l&quot;)</pre>

<p><img alt="plot of chunk unnamed-chunk-14" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-14-1.png?w=604" data-recalc-dims="1"></p>
<h3>Attempt Something that Ultimately Does Not Work</h3>
<p>At first I thought I could just throw two normal distributions at this and it would work. However, it does not work well at all and this kind of not-so-helpful fit turns out to happen a fair bit. I use the <code>nls()</code> function here which is easy to call. It requires a functional form which I set to be:</p>
<p>\[y = C_1 \exp\Big(-{\frac{(x-\mu_1)^2}{2\sigma_1^2}}\Big) + C_2 \exp \Big({-\frac{(x-\mu_2)^2}{2\sigma_2^2}}\Big)\]</p>
<p>where \(\mu_1\) is the \(x\) location of the first peak in \(\gamma\) and \(\mu_2\) is the \(x\) location of the second peak in \(\gamma\). The estimates of \(\sigma_1\) and \(\sigma_2\) can be obtained by trying to estimate the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Full_width_at_half_maximum" target="_blank">full-width-half-maximum</a> (FWHM) of the peaks, which is related to \(\sigma\) by </p>
<p>\[FWHM_i = 2 \sqrt{2\ln2} \times \sigma_i = 2.355 \times \sigma_i\]</p>
<p>I had to first make a little function that returns the respective half-widths at half-maximum and then uses them to estimate the \(FWHM\). Because the peaks are poorly resolved, it also tries to get the smallest possible estimate returning this as <code>FWHM2</code>.</p>

<pre>FWHM.finder &lt;- function(ep.data, mu.index){
  peak.height &lt;- ep.data$y[mu.index]
  fxn.for.roots &lt;- ep.data$y - peak.height/2
  indices &lt;- 1:nrow(ep.data)
  root.indices &lt;- which(diff(sign(fxn.for.roots))!=0)
  tmp &lt;- c(root.indices,mu.index) %>% sort
  tmp2 &lt;- which(tmp == mu.index)
  first.root &lt;- root.indices[tmp2 -1]
  second.root &lt;- root.indices[tmp2]
  HWHM1 &lt;- ep.data$x[mu.index] - ep.data$x[first.root]
  HWHM2 &lt;- ep.data$x[second.root] - ep.data$x[mu.index]
  FWHM &lt;- HWHM2 + HWHM1
  FWHM2 = 2*min(c(HWHM1,HWHM2))
  return(list(HWHM1 = HWHM1,HWHM2 = HWHM2,FWHM = FWHM,FWHM2 = FWHM2))
}</pre>

<p>The peak in the \(\gamma\) region was obtained previously:</p>

<pre>plot(y ~ x, gamma.data, type  = &quot;l&quot;)
gamma.max &lt;- findPeaks(gamma.data$y)
abline(v = gamma.data$x[gamma.max])</pre>

<p><img alt="plot of chunk unnamed-chunk-16" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-16-1.png?w=604" data-recalc-dims="1"></p>
<p>and from them \(\mu_1\) is determined to be 0.7. We have to guess where the second peak is, which is at about \(x=0.75\) and has an index of 252 in the <code>gamma.data</code> dataframe.</p>

<pre>gamma.data[252,]</pre>

<pre>##             x         y
## 252 0.7487757 0.6381026</pre>

<pre>#append the second peak
gamma.max &lt;- c(gamma.max,252)
gamma.mu &lt;- gamma.data$x[gamma.max]
gamma.mu</pre>

<pre>## [1] 0.6983350 0.7487757</pre>

<pre>plot(y ~ x, gamma.data, type  = &quot;l&quot;)
abline(v = gamma.data$x[gamma.max])</pre>

<p><img alt="plot of chunk unnamed-chunk-17" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-17-1.png?w=604" data-recalc-dims="1"></p>
<p>Now we can find the estimates of the standard deviations:</p>

<pre>#find the FWHM estimates of sigma_1 and sigma_2:
FWHM &lt;- lapply(gamma.max, FWHM.finder, ep.data = gamma.data)
gamma.sigma &lt;- unlist(sapply(FWHM, '[', 'FWHM2'))/2.355</pre>

<p>The estimates of \(\sigma_1\) and \(\sigma_2\) are now obtained. The estimates of \(C_1\) and \(C_2\) are just the peak heights.</p>

<pre>peak.heights &lt;- gamma.data$y[gamma.max]</pre>

<p>We can now use <code>nls()</code> to determine the fit.</p>

<pre>fit &lt;- nls(y ~ (C1*exp(-(x-mean1)**2/(2 * sigma1**2)) +
                  C2*exp(-(x-mean2)**2/(2 * sigma2**2))),
           data = gamma.data,
           start = list(mean1 = gamma.mu[1],
                        mean2 = gamma.mu[2],
                        sigma1 = gamma.sigma[1],
                        sigma2 = gamma.sigma[2],
                        C1 = peak.heights[1],
                        C2 = peak.heights[2]),
           algorithm = &quot;port&quot;)</pre>

<p>Determining the fitted values of our unknown coefficients:</p>

<pre>dffit &lt;- data.frame(x=seq(0, 1 , 0.001))
dffit$y &lt;- predict(fit, newdata=dffit)

fit.sum &lt;- summary(fit)
fit.sum #show the fitted coefficients</pre>

<pre>## 
## Formula: y ~ (C1 * exp(-(x - mean1)^2/(2 * sigma1^2)) + C2 * exp(-(x - 
##     mean2)^2/(2 * sigma2^2)))
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## mean1  0.7094793  0.0003312 2142.23   &lt;2e-16 ***
## mean2  0.7813900  0.0007213 1083.24   &lt;2e-16 ***
## sigma1 0.0731113  0.0002382  306.94   &lt;2e-16 ***
## sigma2 0.0250850  0.0011115   22.57   &lt;2e-16 ***
## C1     0.6983921  0.0018462  378.29   &lt;2e-16 ***
## C2     0.0819704  0.0032625   25.12   &lt;2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01291 on 611 degrees of freedom
## 
## Algorithm &quot;port&quot;, convergence message: both X-convergence and relative convergence (5)</pre>

<pre>coef.fit &lt;- fit.sum$coefficients[,1]
mu.fit &lt;- coef.fit[1:2]
sigma.fit &lt;- coef.fit[3:4]
C.fit &lt;- coef.fit[5:6]</pre>

<p>And now we can plot the fitted results against the original results:</p>

<pre>#original
plot(y ~ x, data = gamma.data, type = &quot;l&quot;, main = &quot;This is Garbage&quot;) 
#overall fit
lines(y ~ x, data = dffit, col =&quot;red&quot;, cex = 0.2) 
legend(&quot;topright&quot;, lty = c(1,1,1), col = c(&quot;black&quot;, &quot;green&quot;, &quot;blue&quot;,&quot;red&quot;), c(&quot;Scan&quot;, &quot;Monoclonal&quot;, &quot;Gamma&quot;, &quot;Sum&quot;))
#components of the fit
for(i in 1:2){
  x &lt;- dffit$x
  y &lt;- C.fit[i] *exp(-(x-mu.fit[i])**2/(2 * sigma.fit[i]**2))
  lines(x,y, col = i + 2)
}</pre>

<p><img alt="plot of chunk unnamed-chunk-22" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-22-1.png?w=604" data-recalc-dims="1"></p>
<p><strong>And this is garbage</strong>. The green curve is supposed to be the monoclonal peak, the blue curve is supposed to be the \(\gamma\) background, and the red curve is their sum, the overall fit. This is a horrible failure. </p>
<p>Subsequently, I tried fixing the locations of \(\mu_1\) and \(\mu_2\) but this also yielded similar nonsensical fitting. So, with a lot of messing around trying different functions like the lognormal distribution, the <a rel="nofollow" href="http://pubs.acs.org/doi/abs/10.1021/ac60315a005" target="_blank">Bi-Gaussian</a> distribution and the <a rel="nofollow" href="http://pubs.acs.org/doi/abs/10.1021/ac60319a011?journalCode=ancham" target="_blank">Exponentially Modified Gaussian</a> distribution, and applying various arbitrary weighting functions, and simultaneously fitting the other regions of the SPEP, I concluded that nothing could predictably produce results that represented the clinical reality.</p>
<p>I thought maybe the challenge to obtain a reasonable fit related to the sloping baseline, so I though I would try to remove it. I will model the baseline in the most simplistic manner possible: as a sloped line. </p>
<h2>Baseline Removal</h2>
<p>I will arbitrarily define the tail of the \(\gamma\) region to be those values having \(y \leq 0.02\). Then I will connect the first (x,y) point from the \(\gamma\) region and connect it to the tail. </p>

<pre>gamma.tail &lt;- filter(gamma.data, y &lt;= 0.02) 
baseline.data &lt;- rbind(gamma.data[1,],gamma.tail)
names(baseline.data) &lt;- c(&quot;x&quot;,&quot;y&quot;)
baseline.fun &lt;- approxfun(baseline.data)
plot(y~x, data = gamma.data, type = &quot;l&quot;)
lines(baseline.data$x,baseline.fun(baseline.data$x), col = &quot;blue&quot;)</pre>

<p><img alt="plot of chunk unnamed-chunk-24" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-24-1.png?w=604" data-recalc-dims="1"></p>
<p>Now we can define a new dataframe <code>gamma.no.base</code> that has the baseline removed:</p>

<pre>gamma.no.base &lt;- data.frame(x = gamma.data$x, y = gamma.data$y - baseline.fun(gamma.data$x))
plot(y~x, data = gamma.data, type = &quot;l&quot;)
lines(y ~ x, data = gamma.no.base, lty = 2)
gamma.max &lt;- findPeaks(gamma.no.base$y)[1:2] #rejects a number of extraneous peaks
abline(v = gamma.no.base$x[gamma.max])</pre>

<p><img alt="plot of chunk unnamed-chunk-25" src="https://i0.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-25-1.png?w=604" data-recalc-dims="1"></p>
<p>The black is the original \(\gamma\) and the dashed has the baseline removed. This becomes and easy fit.</p>

<pre>#Estimate the Ci
peak.heights &lt;- gamma.no.base$y[gamma.max]
#Estimate the mu_i
gamma.mu &lt;- gamma.no.base$x[gamma.max] #the same values as before
#Estimate the sigma_i from the FWHM
FWHM &lt;- lapply(gamma.max, FWHM.finder, ep.data = gamma.no.base)
gamma.sigma &lt;- unlist(sapply(FWHM, '[', 'FWHM2'))/2.355

#Perform the fit
fit &lt;- nls(y ~ (C1*exp(-(x-mean1)**2/(2 * sigma1**2)) +
                  C2*exp(-(x-mean2)**2/(2 * sigma2**2))),
           data = gamma.no.base,
           start = list(mean1 = gamma.mu[1],
                        mean2 = gamma.mu[2],
                        sigma1 = gamma.sigma[1],
                        sigma2 = gamma.sigma[2],
                        C1 = peak.heights[1],
                        C2 = peak.heights[2]),
           algorithm = &quot;port&quot;)

#Plot the fit
dffit &lt;- data.frame(x=seq(0, 1 , 0.001))
dffit$y &lt;- predict(fit, newdata=dffit)
fit.sum &lt;- summary(fit)
coef.fit &lt;- fit.sum$coefficients[,1]
mu.fit &lt;- coef.fit[1:2]
sigma.fit &lt;- coef.fit[3:4]
C.fit &lt;- coef.fit[5:6]

plot(y ~ x, data = gamma.no.base, type = &quot;l&quot;)
legend(&quot;topright&quot;, lty = c(1,1,1), col = c(&quot;black&quot;, &quot;green&quot;, &quot;blue&quot;,&quot;red&quot;), c(&quot;Scan&quot;, &quot;Monoclonal&quot;, &quot;Gamma&quot;, &quot;Sum&quot;))
lines(y ~ x, data = dffit, col =&quot;red&quot;, cex = 0.2)
for(i in 1:2){
  x &lt;- dffit$x
  y &lt;- C.fit[i] *exp(-(x-mu.fit[i])**2/(2 * sigma.fit[i]**2))
  lines(x,y, col = i + 2)
}</pre>

<p><img alt="plot of chunk unnamed-chunk-26" src="https://i1.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-26-1.png?w=604" data-recalc-dims="1"></p>
<p>Lo and behold…something that is not completely insane. The green is the monoclonal, the blue is the \(\gamma\) background and the red is their sum, that is, the overall fit. A better fit could now we sought with weighting or with a more flexible distribution shape. In any case, the green peak is now easily determined. Since </p>
<p>\[\int_{-\infty}^{\infty} C_1 \exp\Big(-{\frac{(x-\mu_1)^2}{2\sigma_1^2}}\Big)dx = \sqrt{2\pi}\sigma C_1\]</p>

<pre>A.mono &lt;- sqrt(2*pi)*sigma.fit[1]*C.fit[1] %>% unname() 
A.mono &lt;- round(A.mono,3)
A.mono</pre>

<pre>## sigma1 
##  0.024</pre>

<p>So this peak is 2.4% of the total area. Now, of course, this assumes that nothing under the baseline is attributable to the monoclonal peak and all belongs to normal \(\gamma\)-globulins, which is very unlikely to be true. However, the drop and tangent skimming methods also make assumptions about how the area under the curve contributes to the monoclonal protein. The point is to try to do something that will produce consistent results that can be followed over time. Obviously, if you thought there were three peaks in the \(\gamma\)-region, you'd have to set up your model accordingly. </p>
<h3>All about that Base(line)</h3>
<p>There are obviously better ways to model the baseline because this approach of a linear baseline is not going to work in situations where, for example, there is a small monoclonal in fast \(\gamma\) dwarfed by normal \(\gamma\)-globulins. That is, like this:</p>
<p><img alt="plot of chunk unnamed-chunk-28" src="https://i2.wp.com/labrtorian.com/wp-content/uploads/2017/08/unnamed-chunk-28-1.png?w=604" data-recalc-dims="1"></p>
<p>Something curvilinear or piecewise continuous and flexible enough for more circumstances is generally required. </p>
<p>There is also no guarantee that baseline removal, whatever the approach, is going to be a good solution in other circumstances. Given the diversity of monoclonal peak locations, sizes and shapes, I suspect one would need a few different approaches for different circumstances.</p>
<h2>Conclusions</h2>
<ul>
<li>
<p>The data in the PDFs generated by EP software are processed (probably with splining or similar) followed by the stair-stepping seen above. It would be better to work with raw data from the scanner.</p>
<ul>
<li>This is particularly important if you are using <code>nls()</code> because <code>nls()</code> does not play nice with data having no noise (<a rel="nofollow" href="https://stat.ethz.ch/R-manual/R-devel/library/stats/html/nls.html" target="_blank">“Do not use nls on artificial 'zero-residual' data”</a>)</li>
</ul>
</li>
<li>
<p>Integrating monoclonal peaks under the \(\gamma\) baseline (or \(\beta\)) is unlikely to be a one-size-fits all approach and may require application of a number of strategies to get meaningful results.</p>
<ul>
<li>Basline removal might be helpful at times.</li>
</ul>
</li>
<li>
<p>Peak integration will require human adjudication.</p>
</li>
<li>
<p>While most monoclonal peaks show little skewing, better fitting is likely to be obtained with distributions that afford some skewing.</p>
</li>
<li>
<p><a rel="nofollow" href="http://clinchem.aaccjnls.org/content/early/2016/08/11/clinchem.2015.253740" target="_blank">MASSFIX</a> may soon make this entire discussion irrelevant.</p>
</li>
</ul>
<p><strong>Parting Thought</strong></p>
<p><em>On the matter of fitting</em></p>
<p>In bringing many sons and daughters to glory, it was <em>fitting</em> that God, for whom and through whom everything exists, should make the pioneer of their salvation perfect through what he suffered.</p>
<p>Heb 2:10</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/28/non-linear-regression-application-to-monoclonal-peak-integration-in-serum-protein-electrophoresis/"> 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/TYqkAwUn3us" height="1"></a></div>