14
$\begingroup$

I am trying to create a figure which shows the relationship between viral copies and genome coverage (GCC). This is what my data looks like:

Viral load vs GCC

At first, I just plotted a linear regression but my supervisors told me that was incorrect, and to try a sigmoidal curve. So I did this using geom_smooth:

library(scales)
ggplot(scatter_plot_new, aes(x = Copies_per_uL, y = Genome_cov, colour = Virus)) +
    geom_point() +
    scale_x_continuous(trans = log10_trans(), breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x))) +
        geom_smooth(method = "gam", formula = y ~ s(x), se = FALSE, size = 1) +
    theme_bw() +
    theme(legend.position = 'top', legend.text = element_text(size = 10), legend.title = element_text(size = 12), axis.text = element_text(size = 10), axis.title = element_text(size=12), axis.title.y = element_text(margin = margin (r = 10)), axis.title.x = element_text(margin = margin(t = 10))) +
    labs(x = "Virus copies/µL", y = "GCC (%)") +
    scale_y_continuous(breaks=c(25,50,75,100))

Viral load vs GCC - geom_smooth

However, my supervisors say this is incorrect too because the curves make it look like GCC can go over 100%, which it can't.

My question is: what is the best way to show the relationship between virus copies and GCC? I want to make it clear that A) low virus copies = low GCC, and that B) after a certain amount of virus copies the GCC plateaus.

I've researched a lot of different methods - GAM, LOESS, logistic, piecewise - but I don't know how to tell what is the best method for my data.

EDIT: this is the data:

>print(scatter_plot_new)  
Subsample   Virus   Genome_cov  Copies_per_uL
1   S1.1_RRAV   RRAV    100 92500
2   S1.2_RRAV   RRAV    100 95900
3   S1.3_RRAV   RRAV    100 92900
4   S2.1_RRAV   RRAV    100 4049.54
5   S2.2_RRAV   RRAV    96.9935 3809
6   S2.3_RRAV   RRAV    94.5054 3695.06
7   S3.1_RRAV   RRAV    3.7235  86.37
8   S3.2_RRAV   RRAV    11.8186 84.2
9   S3.3_RRAV   RRAV    11.0929 95.2
10  S4.1_RRAV   RRAV    0   2.12
11  S4.2_RRAV   RRAV    5.0799  2.71
12  S4.3_RRAV   RRAV    0   2.39
13  S5.1_RRAV   RRAV    4.9503  0.16
14  S5.2_RRAV   RRAV    0   0.08
15  S5.3_RRAV   RRAV    4.4147  0.08
16  S1.1_UMAV   UMAV    5.7666  1.38
17  S1.2_UMAV   UMAV    26.0379 1.72
18  S1.3_UMAV   UMAV    7.4128  2.52
19  S2.1_UMAV   UMAV    21.172  31.06
20  S2.2_UMAV   UMAV    16.1663 29.87
21  S2.3_UMAV   UMAV    9.121   32.82
22  S3.1_UMAV   UMAV    92.903  627.24
23  S3.2_UMAV   UMAV    83.0314 615.36
24  S3.3_UMAV   UMAV    90.3458 632.67
25  S4.1_UMAV   UMAV    98.6696 11180
26  S4.2_UMAV   UMAV    98.8405 12720
27  S4.3_UMAV   UMAV    98.7939 8680
28  S5.1_UMAV   UMAV    98.6489 318200
29  S5.2_UMAV   UMAV    99.1303 346100
30  S5.3_UMAV   UMAV    98.8767 345100
$\endgroup$
  • 6
    $\begingroup$ Seems like a logistic regression would be best, since this is bounded between 0 and 100%. $\endgroup$ – mkt Jul 22 at 11:03
  • 1
    $\begingroup$ Try (2) piece-wise (linear) model. $\endgroup$ – user158565 Jul 22 at 14:28
  • 3
    $\begingroup$ try adding method.args=list(family=quasibinomial)) in the arguments to geom_smooth() in your original ggplot code. $\endgroup$ – Ben Bolker Jul 22 at 19:41
  • 4
    $\begingroup$ PS I would encourage you not to suppress standard errors with se=FALSE. Always nice to show people how large the uncertainty actually is ... $\endgroup$ – Ben Bolker Jul 22 at 20:12
  • 2
    $\begingroup$ You don't have enough data points in the transition region to claim with any authority that there's a smooth curve. I could just as easily fit a Heaviside function to the points you are showing us. $\endgroup$ – Carl Witthoft Jul 23 at 13:07
6
$\begingroup$

another way to go about this would be to use a Bayesian formulation, it can be a bit heavy going to start with but it tends to make it much easier to express specifics of your problem as well as getting better ideas of where the "uncertainty" is

Stan is a Monte Carlo sampler with a relatively easy to use programmatic interface, libraries are available for R and others but I'm using Python here

we use a sigmoid like everybody else: it has biochemical motivations as well as being mathematically very convenient to work with. a nice parameterization for this task is:

import numpy as np

def sigfn(x, alpha, beta):
    return 1 / (1 + np.exp(-(x - alpha) * beta))

where alpha defines the midpoint of the sigmoid curve (i.e. where it crosses 50%) and beta defines the slope, values nearer zero are flatter

to show what this looks like, we can pull in your data and plot it with:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_table('raw_data.txt', delim_whitespace=True)
df.columns = ['subsample', 'virus', 'coverage', 'copies']
df.coverage /= 100

x = np.logspace(-1, 6, 201)
plt.semilogx(x, sigfn(np.log(x), 5.5, 3), label='sigfn', color='C2')

sns.scatterplot(df.copies, df.coverage, hue=df.virus, edgecolor='none')

where raw_data.txt contains the data you gave and I transformed the coverage to something more useful. the coefficients 5.5 and 3 look nice and give a plot very much like the other answers:

plot data and manual fit

to "fit" this function using Stan we need to define our model using its own language that's a mix between R and C++. a simple model would be something like:

data {
    int<lower=1> N;  // number of rows
    vector[N] log_copies;
    vector<lower=0,upper=1>[N] coverage;
}
parameters {
    real alpha;
    real beta;
    real<lower=0> sigma;
}
model {
    vector[N] mu;
    mu = 1 ./ (1 + exp(-(log_copies - alpha) * beta));

    sigma ~ cauchy(0, 0.1);
    alpha ~ normal(0, 5);
    beta ~ normal(0, 5);

    coverage ~ normal(mu, sigma);
}

which hopefully reads OK. we have a data block that defines the data we expect when we sample the model, parameters define the things that are sampled, and model defines the likelihood function. You tell Stan to "compile" the model, which takes a while, and then you can sample from it with some data. for example:

import pystan

model = pystan.StanModel(model_code=code)
model.sampling(data=dict(
    N=len(df),
    log_copies=np.log(df.copies),
    coverage=df.coverage,
), iter=10000, chains=4, thin=10)

import arviz
arviz.plot_trace(fit)

arviz makes nice diagnostic plots easy, while printing the fit gives you a nice R-style parameter summary:

4 chains, each with iter=10000; warmup=5000; thin=10; 
post-warmup draws per chain=500, total post-warmup draws=2000.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
alpha   5.51  6.0e-3   0.26   4.96   5.36   5.49   5.64   6.12   1849    1.0
beta    2.89    0.04   1.71   1.55   1.98   2.32   2.95   8.08   1698    1.0
sigma   0.08  2.7e-4   0.01   0.06   0.07   0.08   0.09    0.1   1790    1.0
lp__   57.12    0.04   1.76   52.9   56.1  57.58  58.51  59.19   1647    1.0

the large standard deviation on beta says that the data really doesn't provide much information about this parameter. also some of the answers giving 10+ significant digits in their model fits are overstating things somewhat

because some answers noted that each virus might need its own parameters I extended the model to allow alpha and beta to vary by "Virus". it all gets a bit fiddly, but the two viruses almost certainly have different alpha values (i.e. you need more copies/μL of RRAV for the same coverage) and a plot showing this is:

plot of data and MC samples

the data is the same as before, but I've drawn a curve for 40 samples of the posterior. UMAV seems relatively well determined, while RRAV could follow the same slope and need a higher copy count, or have a steeper slope and a similar copy count. most of the posterior mass is on needing a higher copy count, but this uncertainty might explain some of the differences in other answers finding different things

I mostly used answering this as an exercise to improve my knowledge of Stan, and I've put a Jupyter notebook of this here in case anyone is interested/wants to replicate this.

$\endgroup$
13
$\begingroup$

(Edited taking into account comments below. Thanks to @BenBolker & @WeiwenNg for helpful input.)

Fit a fractional logistic regression to the data. It is well suited to percentage data that is bounded between 0 and 100% and is well-justified theoretically in many areas of biology.

Note that you might have to divide all values by 100 to fit it, since programs frequently expect the data to range between 0 and 1. And as Ben Bolker recommends, to address possible problems caused by the binomial distribution's strict assumptions regarding variance, use a quasibinomial distribution instead.

I've made some assumptions based on your code, such as that there are 2 viruses you are interested in and they may show different patterns (i.e. there may be an interaction between virus type and number of copies).

First, the model fit:

dat <- read.csv('Book1.csv')
dat$logcopies <- log10(dat$Copies_per_uL)
dat$Genome_cov_norm <- dat$Genome_cov/100

fit <- glm(Genome_cov_norm ~ logcopies * Virus, data = dat, family = quasibinomial())
summary(fit)


Call:
glm(formula = Genome_cov_norm ~ logcopies * Virus, family = quasibinomial(), 
    data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.55073  -0.13362   0.07825   0.20362   0.70086  

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)  
(Intercept)          -5.9702     2.8857  -2.069   0.0486 *
logcopies             2.3262     1.0961   2.122   0.0435 *
VirusUMAV             2.6147     3.3049   0.791   0.4360  
logcopies:VirusUMAV  -0.6028     1.3173  -0.458   0.6510  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.6934319)

    Null deviance: 30.4473  on 29  degrees of freedom
Residual deviance:  2.7033  on 26  degrees of freedom

If you trust the p-values, the output does not suggest that the two viruses differ meaningfully. This is in contrast to @NickCox's results below, though we used different methods. I'd not be very confident either way with 30 data points.

Second, the plotting:

It's not hard to code up a way to visualize the output yourself, but there appears to be a ggPredict package that will do most of the work for you (can't vouch for it, I haven't tried it myself). The code will look something like:

library(ggiraphExtra)
ggPredict(fit) + theme_bw(base_size = 20) + geom_line(size = 2) 

Update: I no longer recommend the code or the ggPredict function more generally. After trying it out I found that the plotted points don't exactly reflect the input data but instead are changed for some bizarre reason (some of the plotted points were above 1 and below 0). So I recommend coding it up yourself, though that is more work.

$\endgroup$
  • 6
    $\begingroup$ I endorse this answer, but I'd like to make a point of clarification: I'd call this fractional logistic regression. I think this term would be more widely recognized. When most people hear "logistic regression", I bet they think of a 0/1 dependent variable. One good Stackexchange answer dealing with this nomenclature is here: stats.stackexchange.com/questions/216122/… $\endgroup$ – Weiwen Ng Jul 22 at 15:54
  • 2
    $\begingroup$ @teaelleceecee You evidently have to divide the coverage by 100 first. $\endgroup$ – Nick Cox Jul 22 at 17:49
  • 3
    $\begingroup$ use family=quasibinomial() to avoid the warning (and the underlying problems with too-strict variance assumptions). Take @mkt's advice on the other problem. $\endgroup$ – Ben Bolker Jul 22 at 19:39
  • 2
    $\begingroup$ This may work, but I would like to warn folks that you should have a premise prior to fitting a function that your data in fact should follow that function. Otherwise you are pretty much shooting at random when you pick a fitting function, and you may be fooled by the results. $\endgroup$ – Carl Witthoft Jul 23 at 13:10
  • 5
    $\begingroup$ @CarlWitthoft We hear the sermon but are sinners outside the service. What prior premise led you to suggest a Heaviside function in other comments? The biology here doesn't resemble transition at a sharp threshold. The fact of research here as I understand is that formal theory is weaker than the data. I agree: if people think that a step function makes sense, they should fit one. $\endgroup$ – Nick Cox Jul 23 at 13:20
11
$\begingroup$

This isn't a different answer from @mkt but graphs in particular won't fit into a comment. I first fit a logistic curve in Stata (after logging the predictor) to all data and get this graph

enter image description here

An equation is

100 invlogit(-4.192654 + 1.880951 log10(Copies))

Now I fit curves separately for each virus in the simplest scenario of virus defining an indicator variable. Here for the record is a Stata script:

clear 
input id str9 Subsample   str4 Virus   Genome_cov  Copies_per_uL
1   S1.1_RRAV   RRAV    100 92500
2   S1.2_RRAV   RRAV    100 95900
3   S1.3_RRAV   RRAV    100 92900
4   S2.1_RRAV   RRAV    100 4049.54
5   S2.2_RRAV   RRAV    96.9935 3809
6   S2.3_RRAV   RRAV    94.5054 3695.06
7   S3.1_RRAV   RRAV    3.7235  86.37
8   S3.2_RRAV   RRAV    11.8186 84.2
9   S3.3_RRAV   RRAV    11.0929 95.2
10  S4.1_RRAV   RRAV    0   2.12
11  S4.2_RRAV   RRAV    5.0799  2.71
12  S4.3_RRAV   RRAV    0   2.39
13  S5.1_RRAV   RRAV    4.9503  0.16
14  S5.2_RRAV   RRAV    0   0.08
15  S5.3_RRAV   RRAV    4.4147  0.08
16  S1.1_UMAV   UMAV    5.7666  1.38
17  S1.2_UMAV   UMAV    26.0379 1.72
18  S1.3_UMAV   UMAV    7.4128  2.52
19  S2.1_UMAV   UMAV    21.172  31.06
20  S2.2_UMAV   UMAV    16.1663 29.87
21  S2.3_UMAV   UMAV    9.121   32.82
22  S3.1_UMAV   UMAV    92.903  627.24
23  S3.2_UMAV   UMAV    83.0314 615.36
24  S3.3_UMAV   UMAV    90.3458 632.67
25  S4.1_UMAV   UMAV    98.6696 11180
26  S4.2_UMAV   UMAV    98.8405 12720
27  S4.3_UMAV   UMAV    98.7939 8680
28  S5.1_UMAV   UMAV    98.6489 318200
29  S5.2_UMAV   UMAV    99.1303 346100
30  S5.3_UMAV   UMAV    98.8767 345100
end 

gen log10Copies = log10(Copies)
gen Genome_cov_pr = Genome_cov / 100
encode Virus, gen(virus)
set seed 2803 
fracreg logit Genome_cov_pr log10Copies i.virus, vce(bootstrap, reps(10000)) 

twoway function invlogit(-5.055519 + 1.961538 * x), lc(orange) ra(log10Copies)      ///
|| function invlogit(-5.055519 + 1.233273 + 1.961538 * x), ra(log10Copies) lc(blue) ///
|| scatter Genome_cov_pr log10Copies if Virus == "RRAV", mc(orange) ms(Oh)          ///
|| scatter Genome_cov_pr log10Copies if Virus == "UMAV", mc(blue) ms(+)             ///
legend(order(4 "UMAV" 3 "RRAV") pos(11) col(1) ring(0))                             ///
xla(-1 "0.1" 0 "1" 1 "10" 2 "100" 3 "10{sup:3}" 4 "10{sup:4}" 5 "10{sup:5}")        ///
yla(0 .25 "25" .5 "50" .75 "75" 1 "100", ang(h))                                    ///
ytitle(Genome coverage (%)) xtitle(Genome copies / {&mu}L) scheme(s1color) 

This is pushing hard on a tiny dataset but the P-value for virus looks supportive of fitting two curves jointly.

Fractional logistic regression                  Number of obs     =         30
                                                Replications      =     10,000
                                                Wald chi2(2)      =      48.14
                                                Prob > chi2       =     0.0000
Log pseudolikelihood = -6.9603063               Pseudo R2         =     0.6646

-------------------------------------------------------------------------------
              |   Observed   Bootstrap                         Normal-based
Genome_cov_pr |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
  log10Copies |   1.961538   .2893965     6.78   0.000     1.394331    2.528745
              |
        virus |
        UMAV  |   1.233273   .5557609     2.22   0.026     .1440018    2.322544
        _cons |  -5.055519   .8971009    -5.64   0.000    -6.813805   -3.297234
-------------------------------------------------------------------------------

enter image description here

$\endgroup$
3
$\begingroup$

Try sigmoid function. There are many formulations of this shape including a logistic curve. Hyperbolic tangent is another popular choice.

Given the plots, I can't rule out a simple step function either. I'm afraid you will not be able to differentiate between a step function and any number of sigmoid specifications. You don't have any observations where your percentage is in 50% range, so the simple step formulation can be the most parsimoinous choice that performs no worse than more complex models

$\endgroup$
  • $\begingroup$ It's worth noting the hyperbolic tangent is closely related to the sigmoid function, viz. $\sigma(x)=\frac12\left(1+\tanh\frac{x}{2}\right)$. $\endgroup$ – J.G. Jul 22 at 20:44
  • 2
    $\begingroup$ @J.G. "sigmoid" is a generic term for an S-curve, as far as I'm concerned, but you're right to point to a link between two specifications of a sigmoid $\endgroup$ – Aksakal Jul 22 at 20:49
2
$\begingroup$

Here are the 4PL (4 parameter logistic) fits, both constrained and unconstrained, with the equation as per C.A. Holstein, M. Griffin, J. Hong, P.D. Sampson, “Statistical Method for Determining and Comparing Limits of Detection of Bioassays”, Anal. Chem. 87 (2015) 9795-9801. The 4PL equation is shown in both figures and the parameter meanings are as follows: a = lower asymptote, b = slope factor, c = inflection point, and d = upper asymptote.

Figure 1 constrains a to equal 0% and d to equal 100%:

Fig. 1 Constrained a & d

Figure 2 has no constraints on the 4 parameters in the 4PL equation:

Fig. 2 No constraints

This was fun, I make no pretence of knowing anything biological and it will be interesting to see how it all settles out!

$\endgroup$
  • $\begingroup$ Thank you, this is really helpful. Just wondering, did you do this in MATLAB with the fit function? $\endgroup$ – teaelleceecee Jul 22 at 20:02
  • 1
    $\begingroup$ I used Igor Pro with the user-defined user function shown in the figures. I have used Igor Pro and its predecessor (Igor) since 1988, but lots of other programs can do the curve fitting, e.g., Origin Pro and the very inexpensive Kaleidagraph. And it appears that you have R and (possibly?) access to Matlab, neither of which I know anything about except that they are extremely capable. Best of success with this and I hope you get good news the next time you discuss things with the supervisors! Also, thanks for posting the data! $\endgroup$ – Ed V Jul 22 at 20:11
2
$\begingroup$

I extracted the data from your scatterplot, and my equation search turned up a 3-parameter logistic type equation as a good candidate: "y = a / (1.0 + b * exp(-1.0 * c * x))", where "x" is the log base 10 per your plot. The fitted parameters were a = 9.0005947126706630E+01, b = 1.2831794858584102E+07, and c = 6.6483431489473155E+00 for my extracted data, a fit of the (log 10 x) original data should yield similar results if you re-fit the original data using my values as initial parameter estimates. My parameter values are yielding R-squared = 0.983 and RMSE = 5.625 on the extracted data.

plot

EDIT: Now that the question has been edited to include the actual data, here is a plot using the above 3-parameter equation and initial parameter estimates.

plot2

$\endgroup$
  • $\begingroup$ There seems to have been an error in your data extraction: you have a bunch of negative percentage values. Also, your maximum values are at about 90% instead of 100% as in the original plot. You may have everything offset by about 10% for some reason. $\endgroup$ – mkt Jul 22 at 16:27
  • $\begingroup$ Meh - this is semi-manually extracted data, the original data is required. This is usually sufficient for equation searches, and of course not for final results - which is why I said to use my extract-o-fit parameter values as initial parameter estimates on the original data. $\endgroup$ – James Phillips Jul 22 at 16:32
  • $\begingroup$ Please note that as the actual data has now been added to the post, I have updated this answer using the updated data. $\endgroup$ – James Phillips Jul 22 at 20:30
  • $\begingroup$ Just to reiterate: application of, e.g., a Heaviside function, may yield similar error values. $\endgroup$ – Carl Witthoft Jul 23 at 13:11
  • 1
    $\begingroup$ @JamesPhillips I will attempt to do so (Heaviside --> errorbars or equivalent) $\endgroup$ – Carl Witthoft Jul 23 at 18:39
2
$\begingroup$

Since I had to open my big mouth about Heaviside, here's the results. I set the transition point to log10(viruscopies) = 2.5 . Then I calculated the standard deviations of the two halves of the data set -- that is, the Heaviside is assuming the data on either side has all derivatives = 0 .

RH side std dev = 4.76
LH side std dev = 7.72

Since it turns out there's 15 samples in each batch, the overall std dev is the mean, or 6.24 .

Assuming the "RMSE" quoted in other answers is "RMS error" overall, the Heaviside function would appear to do at least as well as, if not better than, most of the "Z-curve" (borrowed from photographic response nomenclature) fits here.

edit

Useless graph, but requested in comments:

Heaviside curve fit

$\endgroup$
  • $\begingroup$ Woukd you please post a model and scatterplot similarly to what was done in the other answers? I am most curious to see these results and compare. Please also add RMSE and R-squared values for comparison. I personally have never use the Heaviside function and find this very interesting. $\endgroup$ – James Phillips Jul 23 at 23:44
  • $\begingroup$ @JamesPhillips There really isn't anything to show - obviously the scatterplot is the same; all I did was manually select the transition point and take the raw mean of each set of points (left-hand and right-hand). I'm not sure $R^2$ has much meaning here. $\endgroup$ – Carl Witthoft Jul 24 at 13:24
  • $\begingroup$ My meaning was to make a plot similar to those made in the other answers, for the purpose of direct comparison to those answers. $\endgroup$ – James Phillips Jul 24 at 16:27
  • 2
    $\begingroup$ @JamesPhillips you have two wishes left. Choose wisely :-) $\endgroup$ – Carl Witthoft Jul 24 at 18:11
  • $\begingroup$ Thank you for the plot. I observe that in all plots in other answers, the plotted equation follows the curved shape of the data at the top right - your does not, as is the nature of the Heaviside function. This visually appears to contradict your assertion that the Heaviside function would do as well as the equations posted in the other answers - which is why I had previously requested the RMSE and R-squared values, I suspected that the Heaviside function would not follow the shape of the data in this region and might yield worse values for those fit statistics. $\endgroup$ – James Phillips Jul 25 at 15:41

protected by gung Jul 25 at 16:02

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.