What kind of curve (or model) should I fit to my percentage data?

Question

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:

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))

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

Answer 1

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:

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:

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.

Answer 2

(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.

Answer 3

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

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
-------------------------------------------------------------------------------

Answer 4

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

Answer 5

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%:

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

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

Answer 6

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.

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.

Answer 7

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: