ggplot을 사용하여 S 자 곡선을 점에 맞추기

Question

다양한 용량에서 약물 치료의 반응 측정에 대한 간단한 데이터 프레임이 있습니다.

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)

다음 코드를 사용하여 이러한 점을 플로팅하고 원하는 시각화의 기본 기반을 얻을 수 있습니다.

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

이 데이터로 다음으로하고 싶은 작업은 각 약물에 대해 표시된 점에 맞는 S 자 곡선을 플롯에 추가하는 것입니다. 그런 다음이 곡선에 대한 EC50을 계산하고 싶습니다. 내 데이터에 S 자 곡선의 전체 범위가 없을 수도 있다는 것을 알고 있지만 내가 가지고있는 것으로 가능한 최상의 추정치를 얻고 싶습니다. 또한 drug_1의 최종 지점은 S 자 곡선의 예상 추세를 따르지 않지만 약물이 포함 된 용액이 고농도에서 반응을 억제 할 수 있으므로 실제로 예상치 못한 일이 아닙니다 (각 약물이 다른 용액에 있음). 이 점을 데이터에서 제외하고 싶습니다.

내 데이터에 S 자 곡선을 맞추는 단계에서 멈춰 있습니다. 시그 모이 드 곡선을 데이터에 맞추는 다른 솔루션을 살펴 보았지만 작동하지 않는 것 같습니다.

내 문제에 매우 가까운 한 게시물은 다음과 같습니다. (sigmoid) curve fitting glm in r

그것을 바탕으로 나는 시도했다.

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

이로 인해 다음과 같은 오류가 발생하며 기본적으로 직선을 그리는 것으로 보입니다.

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

나는 또한이 링크에서 해결책을 시도했습니다. 이 oxy-Hb 데이터에 S 자 곡선 맞추기

이 경우 다음과 같은 오류가 발생합니다.

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

플롯에 선이 추가되지 않습니다.

이 두 오류를 모두 찾아 보았지만 내 데이터에 맞는 이유를 찾을 수없는 것 같습니다.

어떤 도움이라도 대단히 감사하겠습니다!

Answer 1

내가 코멘트에서 말했듯이, 나는 geom_smooth()매우 쉬운 문제 에만 사용 합니다. 문제가 발생하자마자 nls대신 사용 합니다.

내 대답은 @Duck과 매우 유사하지만 다음과 같은 차이점이 있습니다.

• 비가 중 및 (역 분산) ​​가중치 적용을 모두 보여줍니다.
• 가중치 적용이 작동하도록하려면 nls2약간 더 강력한 알고리즘을 제공하는 패키지 를 사용해야했습니다.
• 내가 사용하는 SSlogis()초기 파라미터 선택을 자동 (자동 시작)를 얻을 수
• 외부에서 모든 예측을 수행 ggplot2한 다음geom_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)
)

EC50이 xmid매개 변수 와 동일하다고 생각합니다 ... 가중 추정치와 비가 중 추정치 사이에 큰 차이가 있습니다.

Answer 2

나는 당신이 원하는 것에 가까운 다음 접근법을 제안 할 것입니다. 나는 또한 binomial가족을 사용하여 데이터에 대한 설정을 시도했지만 0과 1 사이의 값에 대한 몇 가지 문제가 있습니다.이 경우 각각의 비율을 결정하기 위해 추가 변수가 필요합니다. 다음 줄의 코드는 출력을 스케치하기 위해 비선형 근사를 사용합니다.

처음에는 데이터 :

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

비선형 최소 제곱에서 이상적인 매개 변수를 검색하려면 초기 값을 정의해야합니다. 기본 함수와 함께 다음 코드 nls()를 사용하여 초기 값을 얻습니다.

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

이 매개 변수의 초기 접근 방식으로를 사용하여 플롯을 스케치합니다 geom_smooth(). nls()올바른 매개 변수를 찾기 위해 다시 사용 합니다.

#Plot
ggplot(data=df, aes(y=mean_response, x= conc, color = drug)) +
  geom_pointrange(aes(ymax = (mean_response + std_dev), ymin = (mean_response - std_dev))) +
  geom_smooth(data = df[df$drug=='drug_1',],method = "nls", se = FALSE, formula = y ~ a/(1+exp(-b*(x-c))), method.args = list(start = coef(fm1), algorithm='port'), color = "tomato")+ geom_smooth(data = df[df$drug=='drug_2',],method = "nls", se = FALSE,
              formula = y ~ a/(1+exp(-b*(x-c))),
              method.args = list(start = coef(fm0),
                                 algorithm='port'),
              color = "cyan3")

출력 :