Que tipo de curva (ou modelo) devo ajustar aos meus dados de porcentagem?

15

Estou tentando criar uma figura que mostra a relação entre cópias virais e cobertura do genoma (GCC). É assim que meus dados são:

Carga viral vs CCG

No começo, acabei de traçar uma regressão linear, mas meus supervisores disseram que estava incorreto e que tentava uma curva sigmoidal. Então eu fiz isso usando 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))

Carga viral versus GCC - geom_smooth

No entanto, meus supervisores dizem que isso também está incorreto porque as curvas fazem com que o GCC possa ultrapassar 100%, o que não pode.

Minha pergunta é: qual é a melhor maneira de mostrar a relação entre cópias de vírus e o GCC? Quero deixar claro que A) baixa quantidade de vírus = baixo GCC e que B) após uma certa quantidade de vírus copia os platôs do GCC.

Eu pesquisei vários métodos diferentes - GAM, LOESS, logístico, por partes - mas não sei dizer qual é o melhor método para meus dados.

EDIT: estes são os dados:

>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
teaelleceecee
fonte
6
Parece que uma regressão logística seria melhor, pois ela é limitada entre 0 e 100%.
mkt - Restabelece Monica
11
Tente (2) o modelo linear (peça por peça).
user158565
3
tente adicionar method.args=list(family=quasibinomial))os argumentos geom_smooth()no seu código ggplot original.
Ben Bolker
4
PS Eu encorajo você a não suprimir erros padrão com se=FALSE. Sempre bom para mostrar às pessoas o quão grande a incerteza na verdade é ...
Ben Bolker
2
Você não tem pontos de dados suficientes na região de transição para reivindicar com qualquer autoridade que haja uma curva suave. Eu poderia facilmente ajustar uma função Heaviside aos pontos que você está nos mostrando.
Carl Witthoft 23/07/19

Respostas:

6

Outra maneira de fazer isso seria usar uma formulação bayesiana, pode ser um pouco pesada, mas tende a facilitar a expressão de detalhes específicos do seu problema, além de obter melhores idéias sobre onde está a "incerteza". é

Stan é um amostrador de Monte Carlo com uma interface programática relativamente fácil de usar, bibliotecas estão disponíveis para R e outras, mas estou usando Python aqui

usamos um sigmóide como todo mundo: ele tem motivações bioquímicas e é matematicamente muito conveniente de se trabalhar. uma boa parametrização para esta tarefa é:

import numpy as np

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

onde alphadefine o ponto médio da curva sigmóide (ou seja, onde cruza 50%) e betadefine a inclinação, os valores mais próximos de zero são mais planos

para mostrar como é isso, podemos extrair seus dados e plotá-los com:

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

onde raw_data.txtcontém os dados que você forneceu e eu transformei a cobertura em algo mais útil. os coeficientes 5.5 e 3 têm uma boa aparência e apresentam um gráfico muito parecido com as outras respostas:

dados de plotagem e ajuste manual

Para "ajustar" essa função usando Stan, precisamos definir nosso modelo usando sua própria linguagem, que é uma mistura entre R e C ++. um modelo simples seria algo como:

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

que esperançosamente lê OK. temos um databloco que define os dados que esperamos quando parametersamostramos o modelo, definimos as coisas que são amostradas e modeldefine a função de probabilidade. Você diz a Stan para "compilar" o modelo, o que leva um tempo, e então você pode fazer uma amostra dele com alguns dados. por exemplo:

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 facilita as plotagens de diagnóstico, enquanto imprimir o ajuste fornece um resumo resumido dos parâmetros no estilo R:

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

o grande desvio padrão em betadiz que os dados realmente não fornecem muita informação sobre esse parâmetro. Além disso, algumas das respostas que dão mais de 10 dígitos significativos em seus ajustes de modelo estão exagerando um pouco as coisas

porque algumas respostas observaram que cada vírus pode precisar de seus próprios parâmetros, estendi o modelo para permitir alphae betavariar de acordo com o "vírus". tudo fica um pouco complicado, mas os dois vírus quase certamente têm alphavalores diferentes (ou seja, você precisa de mais cópias / μL de RRAV para a mesma cobertura) e um gráfico que mostra isso é:

plotagem de dados e amostras de MC

os dados são os mesmos de antes, mas desenhei uma curva para 40 amostras da parte posterior. UMAVparece relativamente bem determinado, embora RRAVpossa seguir a mesma inclinação e precisar de uma contagem mais alta de cópias ou ter uma inclinação mais acentuada e uma contagem de cópias semelhante. a maior parte da massa posterior está precisando de uma contagem mais alta de cópias, mas essa incerteza pode explicar algumas das diferenças em outras respostas para encontrar coisas diferentes

Eu costumava responder isso como um exercício para melhorar meu conhecimento de Stan, e coloquei um caderno Jupyter disso aqui , caso alguém esteja interessado / queira replicar isso.

Sam Mason
fonte
14

(Editado levando em consideração os comentários abaixo. Agradecemos a @BenBolker & @WeiwenNg pela contribuição útil.)

Ajuste uma regressão logística fracionária aos dados. É adequado para dados percentuais que são delimitados entre 0 e 100% e é justificado teoricamente em muitas áreas da biologia.

Observe que você pode ter que dividir todos os valores por 100 para ajustá-lo, pois os programas freqüentemente esperam que os dados variem entre 0 e 1. E, como recomenda Ben Bolker, para resolver possíveis problemas causados ​​pelas estritas premissas da distribuição binomial em relação à variação, use um distribuição quase-binomial.

Fiz algumas suposições com base no seu código, como a existência de 2 vírus nos quais você está interessado e eles podem mostrar padrões diferentes (por exemplo, pode haver uma interação entre o tipo de vírus e o número de cópias).

Primeiro, o modelo se encaixava:

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

Se você confiar nos valores p, a saída não sugere que os dois vírus diferem significativamente. Isso contrasta com os resultados do @ NickCox abaixo, embora tenhamos usado métodos diferentes. De qualquer maneira, não ficaria muito confiante com 30 pontos de dados.

Segundo, a plotagem:

Não é difícil codificar uma maneira de visualizar a saída você mesmo, mas parece haver um pacote ggPredict que fará a maior parte do trabalho para você (não posso atestar, não tentei por mim mesmo). O código será parecido com:

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

Atualização: não recomendo mais o código ou a função ggPredict de maneira mais geral. Depois de experimentar, descobri que os pontos plotados não refletem exatamente os dados de entrada, mas são alterados por algum motivo bizarro (alguns dos pontos plotados estavam acima de 1 e abaixo de 0). Portanto, recomendo codificá-lo você mesmo, embora isso seja mais trabalhoso.

mkt - Restabelecer Monica
fonte
7
Apoio essa resposta, mas gostaria de esclarecer: eu chamaria isso de regressão logística fracionada. Eu acho que esse termo seria mais amplamente reconhecido. Quando muitas pessoas ouvem "regressão logística", aposto que pensam em uma variável dependente de 0/1. Uma boa resposta Stackexchange lidar com esta nomenclatura está aqui: stats.stackexchange.com/questions/216122/...
Weiwen Ng
2
@teaelleceecee Você evidentemente precisa dividir a cobertura por 100 primeiro.
Nick Cox
4
use family=quasibinomial()para evitar o aviso (e os problemas subjacentes com premissas de variação muito estritas). Siga o conselho de @ mkt sobre o outro problema.
Ben Bolker
2
Isso pode funcionar, mas eu gostaria de avisar às pessoas que você deve ter uma premissa antes de ajustar uma função que seus dados realmente devem seguir essa função. Caso contrário, você estará atirando aleatoriamente ao escolher uma função de ajuste e poderá ser enganado pelos resultados.
Carl Witthoft 23/07/19
6
@CarlWitthoft Ouvimos o sermão, mas somos pecadores fora do serviço. Que premissa anterior o levou a sugerir uma função Heaviside em outros comentários? A biologia aqui não se assemelha à transição em um limiar agudo. O fato de pesquisar aqui, como eu entendo, é que a teoria formal é mais fraca que os dados. Concordo: se as pessoas pensam que uma função de etapa faz sentido, elas devem se encaixar em uma.
23419 Nick Cox
11

Esta não é uma resposta diferente da @mkt, mas os gráficos em particular não cabem em um comentário. Primeiro ajustei uma curva logística no Stata (após registrar o preditor) em todos os dados e obtive esse gráfico

insira a descrição da imagem aqui

Uma equação é

100 invlogit(-4,192654 + 1,880951 log10( Copies))

Agora, ajusto curvas separadamente para cada vírus no cenário mais simples de definição de uma variável indicadora. Aqui para o registro está um script Stata:

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) 

Isso está pressionando fortemente um conjunto de dados minúsculo, mas o valor P para vírus parece favorável ao ajuste de duas curvas em conjunto.

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

insira a descrição da imagem aqui

Nick Cox
fonte
3

Tente a função sigmóide . Existem muitas formulações dessa forma, incluindo uma curva logística. A tangente hiperbólica é outra escolha popular.

Dadas as parcelas, também não posso descartar uma função simples de etapa. Receio que você não consiga diferenciar entre uma função step e qualquer número de especificações sigmóides. Você não tem nenhuma observação em que sua porcentagem está na faixa de 50%; portanto, a formulação de etapas simples pode ser a escolha mais parcimoniosa que não apresenta desempenho pior do que modelos mais complexos

Aksakal
fonte
σ(x)=1 12(1 1+tanhx2)
2
@JG "sigmóide" é um termo genérico para uma curva em S, tanto quanto eu estou preocupado, mas você está certo para apontar para uma ligação entre duas especificações de um sigmóide
Aksakal
2

Aqui estão os ajustes 4PL (4 parâmetros logísticos), restritos e irrestritos, com a equação conforme CA Holstein, M. Griffin, J. Hong, PD Sampson, "Método Estatístico para Determinar e Comparar Limites de Detecção de Bioensaios", Anal . Chem. 87 (2015) 9795-9801. A equação 4PL é mostrada nas duas figuras e os significados dos parâmetros são os seguintes: a = assíntota inferior, b = fator de inclinação, c = ponto de inflexão ed = assíntota superior.

A Figura 1 restringe a para 0% e d para 100%:

Fig. 1 A & d restrito

A Figura 2 não possui restrições nos 4 parâmetros na equação 4PL:

Fig. 2 Sem restrições

Foi divertido, não pretendo saber nada biológico e será interessante ver como tudo se resolve!

Ed V
fonte
Obrigado, isso é realmente útil. Imaginando, você fez isso no MATLAB com a função de ajuste?
teaelleceecee
11
Usei o Igor Pro com a função de usuário definida pelo usuário, mostrada nas figuras. Eu uso o Igor Pro e seu antecessor (Igor) desde 1988, mas muitos outros programas podem fazer o ajuste de curva, por exemplo, Origin Pro e o Kaleidagraph muito barato. E parece que você tem acesso R e (possivelmente?) Ao Matlab, dos quais nada sei, exceto que eles são extremamente capazes. O melhor é o sucesso com isso e espero que você receba boas notícias na próxima vez que discutir as coisas com os supervisores! Além disso, obrigado por postar os dados!
Ed V
2

Extraí os dados do gráfico de dispersão e minha pesquisa de equações gerou uma equação de tipo logístico de 3 parâmetros como um bom candidato: "y = a / (1,0 + b * exp (-1,0 * c * x))", onde " x "é a base de log 10 por seu gráfico. Os parâmetros ajustados foram a = 9.0005947126706630E + 01, b = 1.2831794858584102E + 07 e c = 6.6483431489473155E + 00 para meus dados extraídos, um ajuste dos dados originais (log 10 x) deve produzir resultados semelhantes se você reajustar os dados originais usando meus valores como estimativas de parâmetro inicial. Meus valores de parâmetro estão produzindo R-quadrado = 0,983 e RMSE = 5,625 nos dados extraídos.

enredo

EDIT: Agora que a pergunta foi editada para incluir os dados reais, aqui está um gráfico usando a equação de 3 parâmetros acima e estimativas de parâmetros iniciais.

plot2

James Phillips
fonte
Parece ter havido um erro na extração de dados: você tem vários valores percentuais negativos. Além disso, seus valores máximos estão em cerca de 90%, em vez de 100%, como no gráfico original. Você pode ter tudo compensado em cerca de 10% por algum motivo.
mkt - Reinstala Monica
Meh - trata-se de dados extraídos semi-manualmente, os dados originais são necessários. Isso geralmente é suficiente para pesquisas de equações e, é claro, não para resultados finais - razão pela qual eu disse para usar meus valores de parâmetro extrair o ajuste como estimativas iniciais de parâmetros nos dados originais.
James Phillips
Observe que, como os dados reais foram adicionados à postagem, atualizei esta resposta usando os dados atualizados.
James Phillips
Apenas para reiterar: a aplicação de, por exemplo, uma função Heaviside, pode gerar valores de erro semelhantes.
Carl Witthoft 23/07/19
11
@JamesPhillips Vou tentar fazê-lo (Heaviside -> barras de erro ou equivalente) #
Carl Witthoft
2

Desde que eu tive que abrir minha boca grande sobre Heaviside, aqui estão os resultados. Defino o ponto de transição como log10 (cópias de vírus) = 2,5. Depois, calculei os desvios padrão das duas metades do conjunto de dados - ou seja, o Heaviside está assumindo que os dados de ambos os lados têm todas as derivadas = 0.

Desvio padrão do lado direito = 4.76 Desvio
padrão do lado esquerdo = 7.72

Como acontece que há 15 amostras em cada lote, o desvio padrão geral é a média ou 6,24.

Supondo que o "RMSE" citado em outras respostas seja "erro RMS" no geral, a função Heaviside parece funcionar tanto quanto, se não melhor do que a maioria da "curva Z" (emprestada da nomenclatura da resposta fotográfica) aqui.

editar

Gráfico inútil, mas solicitado nos comentários:

Ajuste da curva Heaviside

Carl Witthoft
fonte
Você poderia postar um modelo e um gráfico de dispersão similar ao que foi feito nas outras respostas? Estou muito curioso para ver esses resultados e comparar. Adicione também valores RMSE e R ao quadrado para comparação. Pessoalmente, nunca usei a função Heaviside e acho isso muito interessante.
James Phillips
@JamesPhillips Realmente não há nada para mostrar - obviamente o gráfico de dispersão é o mesmo; tudo o que fiz foi selecionar manualmente o ponto de transição e obter a média bruta de cada conjunto de pontos (mão esquerda e mão direita). Eu não tenho certeza tem muito significado aqui. R2
Carl Witthoft 24/07/19
Meu significado era fazer um enredo semelhante ao feito nas outras respostas, com o objetivo de comparar diretamente essas respostas.
James Phillips
2
@JamesPhillips você tem mais dois desejos. Escolha sabiamente :-)
Carl Witthoft 24/07/19
Obrigado pela trama. Observo que em todas as plotagens de outras respostas, a equação plotada segue a forma curva dos dados no canto superior direito - a sua não, assim como a natureza da função Heaviside. Visualmente, isso parece contradizer sua afirmação de que a função Heaviside funcionaria tão bem quanto as equações postadas nas outras respostas - e é por isso que eu solicitei anteriormente os valores RMSE e R-quadrado, suspeitei que a função Heaviside não seguisse a forma dos dados nessa região e pode gerar valores piores para as estatísticas de ajuste.
James Phillips