2016-10-10 5 views
1

Mein Modell geht nicht weiter zu den Asymptoten, wenn es in ggplot2 geplottet wird, obwohl es in R-Basisgrafiken funktioniert. In ggplot2, stoppt sie an bestimmten Punkten auf der X-Achse (Bilder unten enthalten), Ich bin 90% sicher, dass das seq() verwandt ist, Daten, die sie am Ende dieses Beitrags geschriebenPlotten drc Modell in ggplot2; Problem mit seq()

ich predict() bin mit zu transformieren Daten aus einem drm (Dosis-Antwort-Paket) Logit-Modell. In den Basisgrafiken sieht die sigmoidale Kurve gut aus, in ggplot2 nicht so sehr.
(Ich glaube nicht, dass Sie es brauchen, aber nur für den Fall, werden die Daten am Ende dieser Beitrag enthalten)

library(drc) 
library(ggplot2) 

Passdaten Logitmodell

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

Extrahieren der Daten des Modells mit Code aus der Demonstration des Autors (unten verlinkt), habe ich:

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] 

Und schließlich die Grafik 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 zeigen Unterschiede in meinem Grundstück auf seq abhängig.

Wenn das folgende für seq() verwendet wird, wird die Grafik über die Daten hinaus gezogen! (seq(log(0.5),**log(10000)**,length=200))) (Bilder 1 & 2)

Ich verstehe nicht seq() trotz es zu erforschen. Würde jemand bitte mir helfen zu verstehen, was mit meinen Plots los ist?

Es scheint, dass der erste Term in seq die untere Grenze definiert. Aber was definiert dann der dritte Begriff? Sie können dies in Bildern sehen 3 & 4 - Die Grafik ist anständig; Ich habe das Thema etwas verschleiert, aber es geht immer noch nicht in Richtung Infin. Das ist ein kleines Problem, weil ich 8 Logit-Modelle zusammenstellen werde.

- (Ich kann nicht mehr als zwei Links veröffentlichen, so kahl mit diesen Titeln)

Für alle Mühe Plotten der Demokratischen Republik Kongo/drm Modelle mit ggplot2 aufweisen, wobei die folgenden Beiträge waren sehr hilfreich: Die Suche für Plotten-Dosis-Wirkungs-Kurve-mit-ggplot2-and-drc
und diesen Titel: Plotten-Dosis-Wirkungs-Kurve-mit-ggplot2-and-drc

ich habe gefolgt der Autor der Demokratischen Republik Kongo s Anweisungen, die in den unterstützenden Informationen seines Artikels gefunden werden können - Teil dieses Codes wurde oben verwendet. Artikeltitel: Dosis-Antwort-Analyse mit R, Christopher Ritz. Plus eins.

Wenn diese Daten in eine weniger als ideale Format ist, lassen Sie mich das Know-

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

Antwort

1

Dies ist ein langer Kommentar. Ich glaube, Sie haben dose Werte verwechselt, um predict() oder aes(x) Werte zu geben.

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

Siehe ?seq für eine Erklärung dessen, was die Bedingungen seq tun. seq(log(.5), log(1000), length=200) macht 200 Zahlen gleichmäßig im Abstand von Log (0,5) bis Log (1000). Wenn Sie das dritte Argument nicht benennen (oder by=XYZ angeben), ist dies der Abstand zwischen den Zahlen.

Also, wenn Sie seq(log(0.5), log(1000), length=200) tun, wird es Ihre Anpassung an 200 Punkte von Protokoll (0,5) zu Protokoll (1000) berechnen.

Ich denke, mit "in die Unendlichkeit hinausgehen" meinen Sie, dass Sie wollen, dass die Linie über den Rand der Handlung hinausgeht, anstatt vor der Kante anzuhalten, wie in Ihren verknüpften Bildern. Standardmäßig versucht ggplot sicherzustellen, dass all Ihre geplotteten Objekte in Ihr Diagramm passen, sodass die Achsen ein wenig über das Ausmaß Ihrer Daten hinaus erweitert werden.

Wenn Sie es einschränken möchten, verwenden Sie einfach + xlim(c(lower, upper)).

(ich Probleme habe drc Installation so kann Ihr Beispiel nicht reproduzieren; hier ist ein Spielzeug ein)

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

In der oben geht die Linie, um 100 (wie erwartet) und die Grenzen gehen ein wenig darüber hinaus. Wenn ich die Linie wollte den Rand des Grundstücks berühren, dann kann ich (zB) bei 100 genau die x Grenzen Clip - mit dem limx Argumente coord_trans:

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

clipped

Also, wenn Sie Grundstück Entscheiden Sie, was die Grenzen Ihrer x-Achse sein werden und stellen Sie sicher, dass Sie alle Ihre Modelle entlang dieser Werte vorhersagen. Beschränke dann das x-Limit auf diese Grenzen und die Linien scheinen "ins Unendliche zu gehen".

+0

danke für die Erklärung, aber ... ich ein anderes Bild hier gehostet haben [Imgur link] (https://imgur.com/r4xmoB2) Die Die ersten beiden (A1, A2) teilen sich die gleichen 'seq()' Werte, wobei der mittlere Term auf log (10000) gesetzt ist. Die unteren beiden (B1, B2) sind auf log (1000) eingestellt. Ich habe identische 'cartesian_coord()' Limits auf A2, B2 gesetzt. Abgesehen davon gibt es keine weiteren Unterschiede im Code. - Aber die Kurven (zwischen A und B) sind so verschieden. ** Können Sie mir helfen zu verstehen, warum die A- und B-Gruppe sich so anders verhalten? ** – Arch

+0

'geom_point()' scheint die Punkte mit einem Datensatz zu zeichnen, der nicht vom 'seq()' -Begriff betroffen ist; Die Modellkurven sind jedoch immer noch mit niedrigeren Dosierungen auf der X-Achse ausgerichtet, obwohl ich log x 'coord_trans()' command implementiert – Arch

+0

'cartesian_coord' ist ein nicht transformiertes Koordinatensystem, dh Ihre Log-Transformation wird rückgängig gemacht, wenn Sie sie anwenden. Deshalb. Verwenden Sie das Argument 'limx' für' coord_trans', um mit dem Code in Ihrer Frage kompatibel zu sein. –