Ajustando uma curva sigmoidal a pontos com ggplot

Aug 25 2020

Tenho um quadro de dados simples para as medições de resposta de um tratamento com drogas em várias doses:

drug <- c("drug_1", "drug_1", "drug_1", "drug_1", "drug_1", 
  "drug_1", "drug_1", "drug_1", "drug_2", "drug_2", "drug_2", 
        "drug_2", "drug_2", "drug_2", "drug_2", "drug_2")

conc <- c(100.00, 33.33, 11.11, 3.70, 1.23, 0.41, 0.14, 
        0.05, 100.00, 33.33, 11.11, 3.70, 1.23, 0.41, 0.14, 0.05)

mean_response <- c(1156, 1833, 1744, 1256, 1244, 1088, 678, 489, 
        2322, 1867, 1333, 944, 567, 356, 200, 177)

std_dev <- c(117, 317, 440, 200, 134, 38, 183, 153, 719,
      218, 185, 117, 166, 167, 88, 50)

df <- data.frame(drug, conc, mean_response, std_dev)

Posso plotar esses pontos usando o código a seguir e obter a base básica da visualização que gostaria:

p <- ggplot(data=df, aes(y=mean_response, x= conc, color = drug)) +
  geom_pointrange(aes(ymax = (mean_response + std_dev), ymin = (mean_response - std_dev))) +
  scale_x_log10()

p

A próxima coisa que gostaria de fazer com esses dados é adicionar uma curva sigmoidal ao gráfico, que se ajusta aos pontos traçados para cada medicamento. Em seguida, gostaria de calcular o EC50 para esta curva. Percebo que posso não ter toda a gama da curva sigmoidal em meus dados, mas espero obter a melhor estimativa possível com o que tenho. Além disso, o ponto final para a droga_1 não segue a tendência esperada de uma curva sigmoidal, mas isso não é realmente inesperado, pois as soluções em que a droga está podem inibir as respostas em altas concentrações (cada droga está em uma solução diferente). Eu gostaria de excluir este ponto dos dados.

Estou ficando preso na etapa de ajustar uma curva sigmoidal aos meus dados. Eu examinei algumas outras soluções para ajustar curvas sigmoidais aos dados, mas nenhuma parece funcionar.

Um post que está muito próximo do meu problema é este: (sigmóide) ajuste de curva glm em r

Com base nisso, tentei:

p + geom_smooth(method = "glm", family = binomial, se = FALSE)

Isso dá o seguinte erro e parece ser o padrão para traçar linhas retas:

`geom_smooth()` using formula 'y ~ x'
Warning message:
Ignoring unknown parameters: family

Eu também tentei a solução deste link: Ajustando uma curva sigmoidal para esses dados de oxi-Hb

Nesse caso, recebo o seguinte erro:

Computation failed in `stat_smooth()`:
Convergence failure: singular convergence (7)

e nenhuma linha é adicionada ao gráfico.

Eu tentei procurar esses dois erros, mas não consigo encontrar um motivo que faça sentido com meus dados.

Qualquer ajuda seria muito apreciada!

Respostas

2 BenBolker Aug 25 2020 at 06:27

Como disse em um comentário, só usaria geom_smooth()para um problema muito fácil; assim que tiver problemas, uso nls.

Minha resposta é muito semelhante à de @Pato, com as seguintes diferenças:

Mostro ajustes não ponderados e ponderados (variação inversa).
Para fazer os ajustes ponderados funcionarem, tive que usar o nls2pacote, que fornece um algoritmo um pouco mais robusto
Eu uso SSlogis()para obter a seleção automática de parâmetro inicial (auto-inicialização)
Eu faço todas as previsões fora de ggplot2, em seguida, coloco emgeom_line()

p1 <- nls(mean_response~SSlogis(conc,Asym,xmid,scal),data=df,
          subset=(drug=="drug_1" & conc<100)
        ## , weights=1/std_dev^2  ## error in qr.default: NA/NaN/Inf ...
          )

library(nls2)
p1B <- nls2(mean_response~SSlogis(conc,Asym,xmid,scal),data=df,
            subset=(drug=="drug_1" & conc<100),
            weights=1/std_dev^2)

p2 <- update(p1,subset=(drug=="drug_2"))
p2B <- update(p1B,subset=(drug=="drug_2"))

pframe0 <- data.frame(conc=10^seq(log10(min(df$conc)),log10(max(df$conc)), length.out=100))
pp <- rbind(
    data.frame(pframe0,mean_response=predict(p1,pframe0),
               drug="drug_1",wts=FALSE),
    data.frame(pframe0,mean_response=predict(p2,pframe0),
               drug="drug_2",wts=FALSE),
    data.frame(pframe0,mean_response=predict(p1B,pframe0),
               drug="drug_1",wts=TRUE),
    data.frame(pframe0,mean_response=predict(p2B,pframe0),
               drug="drug_2",wts=TRUE)
)

library(ggplot2); theme_set(theme_bw())
(ggplot(df,aes(conc,mean_response,colour=drug)) +
 geom_pointrange(aes(ymin=mean_response-std_dev,
                     ymax=mean_response+std_dev)) +
 scale_x_log10() +
 geom_line(data=pp,aes(linetype=wts),size=2)
)

Eu acredito que EC50 é equivalente ao xmidparâmetro ... observe as grandes diferenças entre as estimativas ponderadas e não ponderadas ...

1 Duck Aug 25 2020 at 05:27

Eu sugeriria a próxima abordagem que está próxima do que você deseja. Também tentei com uma configuração para seus dados usando binomialfamília, mas há alguns problemas sobre valores entre 0 e 1. Nesse caso, você precisaria de uma variável adicional para determinar as respectivas proporções. O código nas linhas a seguir usa uma aproximação não linear para esboçar sua saída.

Inicialmente, os dados:

library(ggplot2)
#Data
df <- structure(list(drug = c("drug_1", "drug_1", "drug_1", "drug_1", 
"drug_1", "drug_1", "drug_1", "drug_1", "drug_2", "drug_2", "drug_2", 
"drug_2", "drug_2", "drug_2", "drug_2", "drug_2"), conc = c(100, 
33.33, 11.11, 3.7, 1.23, 0.41, 0.14, 0.05, 100, 33.33, 11.11, 
3.7, 1.23, 0.41, 0.14, 0.05), mean_response = c(1156, 1833, 1744, 
1256, 1244, 1088, 678, 489, 2322, 1867, 1333, 944, 567, 356, 
200, 177), std_dev = c(117, 317, 440, 200, 134, 38, 183, 153, 
719, 218, 185, 117, 166, 167, 88, 50)), class = "data.frame", row.names = c(NA, 
-16L))

Em um mínimo de quadrados não linear, você precisa definir valores iniciais para a busca de parâmetros ideais. Usamos o próximo código com função de base nls()para obter esses valores iniciais:

#Drug 1
fm1 <- nls(log(mean_response) ~ log(a/(1+exp(-b*(conc-c)))), df[df$drug=='drug_1',], start = c(a = 1, b = 1, c = 1)) #Drug 2 fm2 <- nls(log(mean_response) ~ log(a/(1+exp(-b*(conc-c)))), df[df$drug=='drug_2',], start = c(a = 1, b = 1, c = 1))

Com essa abordagem inicial de parâmetros, esboçamos o gráfico usando geom_smooth(). Novamente usamos nls()para encontrar os parâmetros corretos: