2016-10-10 7 views
1

Моя модель не продолжается по отношению к асимптотам при построении графика в ggplot2, хотя она и есть в базовой графике R. В ggplot2, он останавливается в определенных точках на оси X (включая изображения ниже), Я 90% уверен, что это связано с seq(), данные, которые он разместил в нижней части этого постаПланирование модели drc в ggplot2; issue with seq()

Я использую predict() для преобразования данные из логит-модели drm (доза-реакция). В базовой графике сигмоидальная кривая выглядит великолепно, в ggplot2, не так много.
(я не думаю, что вам это нужно, но только в том случае, данные включены в нижней части этого поста)

library(drc) 
library(ggplot2) 

установочные данные в логит модели

mod1 <- drm(probability ~ (dose), weights = total, data = mydata1,  type ="binomial", fct = LL.2()) 

plot(mod1,broken=FALSE,type="all",add=FALSE, col= "purple", xlim=c(0, 10000)) 

Image of Base graphic 2-parameter logit

Извлечение данных модели с использованием кода из демонстрации автора (см. Ниже), у меня есть:

newdata1 <-expand.grid(dose=exp(seq(log(0.5),log(100),length=200))) 

pm1<- predict(mod1, newdata=newdata1,interval="confidence") 

newdata1$p1 <-pm1[,1] 

newdata1$pmin1 <-pm1[,2] 

newdata1$pmax1 <- pm1[,3] 

И, наконец, графика ggplot2.

p1 <- ggplot(mydata1, aes(x=dose01,y=probability))+ 
    geom_point()+ 
    geom_ribbon(data=newdata1, aes(x=dose,y=p1, 
ymin=pmin1,ymax=pmax1),alpha=0.2,color="blue",fill="pink") + 
    geom_step(data=newdata1, aes(x=dose,y=p1))+ 
    coord_trans(x="log") + #creates logline for x axis 
    xlab("dose")+ylab("response") 

Images 1&2 and 3&4 показать отличия в моем участке в зависимости от сл.

Если для seq() использовано следующее, графика будет удалена из-за данных! (seq(log(0.5),**log(10000)**,length=200))) (изображения 1 & 2)

Я не понимаю seq() несмотря на изучение его. Хотел бы кто-нибудь помочь мне понять, что происходит с моими сюжетами?

Похоже, что первый член в seq определяет нижнюю границу. но тогда какой третий термин определяет? Вы можете видеть это на изображениях 3 & 4 - Графика приличная; Я несколько запутал проблему, но она по-прежнему не продолжается до бесконечности. Это небольшая проблема, потому что я буду совмещать 8 моделей логита.

- (я не могу отправить более чем две ссылки, так что голые с этими названиями)

Для тех, кто возникли проблемы построения модели DRC/DrM с ggplot2 следующие сообщения были очень полезны: Поиск для Plotting-доза-реакция-кривые-с-ggplot2-и-дрк
и это название: черчения-доза-реакция-кривые-с-ggplot2-и-дрк

Я следовал автор ДРК " которые могут быть найдены в вспомогательной информации его статьи. Часть этого кода использовалась выше. Название статьи: анализ доза-ответа с использованием R, Кристофер Ритц. PlosOne.

Если данные находятся в менее чем идеальный формат, позвольте мне ноу

> dput(mydata1) 

structure(list(dose = c(25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 75L, 75L, 75L, 75L, 75L, 75L, 
75L, 75L, 75L, 75L, 75L, 75L, 75L, 75L, 75L, 100L, 100L, 100L, 
100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 100L, 
100L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 150L, 
150L, 150L, 150L, 150L, 150L, 200L, 200L, 200L, 200L, 200L, 200L, 
200L, 200L, 200L, 200L, 200L, 200L, 200L, 200L, 200L), total = c(25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 
25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L, 25L), affected = c(2L, 
0L, 0L, 0L, 0L, 1L, 1L, 2L, 2L, 4L, 1L, 1L, 4L, 0L, 10L, 0L, 
1L, 0L, 1L, 0L, 3L, 0L, 4L, 2L, 0L, 2L, 0L, 1L, 2L, 3L, 2L, 0L, 
2L, 0L, 0L, 4L, 0L, 1L, 2L, 3L, 0L, 21L, 1L, 3L, 1L, 2L, 7L, 
0L, 0L, 0L, 0L, 8L, 7L, 3L, 7L, 2L, 2L, 10L, 3L, 4L, 0L, 7L, 
0L, 3L, 3L, 20L, 25L, 22L, 23L, 22L, 18L, 14L, 20L, 20L, 21L), 
    probability = c(0.08, 0, 0, 0, 0, 0.04, 0.04, 0.08, 0.08, 
    0.16, 0.04, 0.04, 0.16, 0, 0.4, 0, 0.04, 0, 0.04, 0, 0.12, 
    0, 0.16, 0.08, 0, 0.08, 0, 0.04, 0.08, 0.12, 0.08, 0, 0.08, 
    0, 0, 0.16, 0, 0.04, 0.08, 0.12, 0, 0.84, 0.04, 0.12, 0.04, 
    0.08, 0.28, 0, 0, 0, 0, 0.32, 0.28, 0.12, 0.28, 0.08, 0.08, 
    0.4, 0.12, 0.16, 0, 0.28, 0, 0.12, 0.12, 0.8, 1, 0.88, 0.92, 
    0.88, 0.72, 0.56, 0.8, 0.8, 0.84)), .Names = c("dose", "total", 
"affected", "probability"), row.names = c(NA, -75L), class = "data.frame") 

ответ

1

Это длинный комментарий. Я думаю, вы ошиблись в значениях dose, чтобы дать predict() или aes(x) значения.

log10000 <- exp(seq(log(0.5), log(10000), length=200)) 
log1000 <- exp(seq(log(0.5), log(1000), length=200)) 

log10000df <- as.data.frame(cbind(dose = log10000, predict(mod1, data.frame(dose = log10000), interval="confidence"))) 
log1000df <- as.data.frame(cbind(dose = log1000, predict(mod1, data.frame(dose = log1000), interval="confidence"))) 

## a common part 
p0 <- ggplot(mydata1, aes(x = dose, y = probability)) + 
    geom_point() + coord_trans(x="log") + 
    xlab("dose") + ylab("response") + xlim(0.5, 10001) 

p10000 <- p0 + geom_line(data = log10000df, aes(x = dose, y = Prediction)) + 
    geom_ribbon(data = log10000df, aes(x = dose, y = Prediction, ymin = Lower, ymax = Upper), 
       alpha = 0.2, color = "blue", fill = "pink") 

p1000 <- p0 + geom_line(data = log1000df, aes(x = dose, y = Prediction)) + 
    geom_ribbon(data = log1000df, aes(x = dose, y = Prediction, ymin = Lower, ymax = Upper), 
       alpha = 0.2, color = "blue", fill = "pink") 

enter image description here enter image description here

0

См ?seq для объяснения того, что условия seq сделать. seq(log(.5), log(1000), length=200) делает 200 номеров равномерно разнесенными, начиная с журнала (0,5), идущего в журнал (1000). Если вы не укажете третий аргумент (или укажите by=XYZ), это будет интервал между числами.

Итак, когда вы сделаете seq(log(0.5), log(1000), length=200), он рассчитает вашу посадку в 200 точках от журнала (0,5) до журнала (1000).

Я думаю, что «выходите на бесконечность» вы подразумеваете, что вы хотите, чтобы линия уходила с края графика, а не останавливалась перед краем, как в ваших связанных изображениях. По умолчанию ggplot будет пытаться обеспечить, чтобы все ваши построенные объекты соответствовали вашему сюжету, поэтому он значительно расширяет оси, не ограничиваясь объемами ваших данных.

Если вы хотите ограничить его, просто используйте + xlim(c(lower, upper)).

(я имею проблемы с установкой drc поэтому не может воспроизвести ваш пример, вот игрушка один)

x = seq(0.5, 100, length=200) 
df <- data.frame(x=x, y=x^2) 
ggplot(df, aes(x=x, y=y)) + geom_line() + coord_trans(x="log") 

original plot

В приведенном выше, строка выходит до 100 (как и ожидалось) и пределы простираются немного дальше. Если бы я хотел линию, чтобы прикоснуться к краю участка, то я могу (например) обрезают пределы х на 100 точно - использовать limx аргумент coord_trans:

ggplot(df, aes(x=x, y=y)) + geom_line() + coord_trans(x="log", limx=c(0.5, 100)) 

clipped

Итак, когда вы сюжет ваши модели, определите, какие будут границы вашей оси x, и убедитесь, что вы прогнозируете все ваши модели по этим значениям. Затем ограничьте предел x этими границами, и линии появятся «перейти на бесконечность».

+0

спасибо за объяснение, но ... я принимал другое изображение здесь [Imgur ссылка] (https://imgur.com/r4xmoB2) верхние два (A1, A2) используют одни и те же значения 'seq()', при этом средний термин устанавливается в log (10000). Нижние два (B1, B2) установлены на log (1000). Я поставил идентичные 'cartesian_coord()' ограничения на A2, B2. Помимо этого, в коде нет других отличий. - Но кривые (между А и В) настолько различны. ** Можете ли вы помочь мне понять, почему группы A и B ведут себя по-разному? ** – Arch

+0

'geom_point()' похоже, выводит точки с использованием набора данных, не зависящего от члена 'seq()'; однако модельные кривые по-прежнему выравниваются с более низкими дозами по оси X, несмотря на то, что я реализую log x 'coord_trans()' команда – Arch

+0

'cartesian_coord' - это нетрансформированная система координат, то есть ваш логарифмический перевод отменяется, когда вы его применяете. Вот почему. Используйте аргумент 'limx' для' coord_trans' для совместимости с кодом в вашем вопросе. –