Polynomiale (polynomielle) Regression mit poly
Moderator: EDi
Re: Polynomiale (polynomielle) Regression mit poly
Hallo Andreas,
mein erstes Beispiel zeigt, dass und wie man von ggplot2 Regressionkurven von polynomialen Regressionen anzeigen lassen kann. Das dient dann aber nur der Anzeige von Kurven und ist m. W. sonst nicht zu gebrauchen (um z. B. Werte vorherzusagen oder Signifikanztests durchzuführen).
Mein zweites Beispiel zeigt, wie man aus mit lm durchgeführte Regressionen Werte vorhersagen kann und wie man damit Kurven zeichnen kann. Das erfordert einen Umweg über lm und predict, dafür erhält man dann auch Werte, die man außerhalb von Abbildungen verwenden kann.
LG,
Bernhard
mein erstes Beispiel zeigt, dass und wie man von ggplot2 Regressionkurven von polynomialen Regressionen anzeigen lassen kann. Das dient dann aber nur der Anzeige von Kurven und ist m. W. sonst nicht zu gebrauchen (um z. B. Werte vorherzusagen oder Signifikanztests durchzuführen).
Mein zweites Beispiel zeigt, wie man aus mit lm durchgeführte Regressionen Werte vorhersagen kann und wie man damit Kurven zeichnen kann. Das erfordert einen Umweg über lm und predict, dafür erhält man dann auch Werte, die man außerhalb von Abbildungen verwenden kann.
LG,
Bernhard
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Re: Polynomiale (polynomielle) Regression mit poly
So, jetzt habe ich mich auch noch mal am Aufhübschen Deiner Plots versucht - Hausarrest macht's möglich !
Ich habe zwei in meinen Augen wesentliche Änderungen drin: zunächst habe ich den Plotdatensatz ins "lange Format" umgewandelt -
Zudem gewinnst Du die Freiheit, bei Bedarf die Grafiken im Handumdrehen über "facets" zu entflechten...
Zum zweiten habe ich die Anmerkungen zu den verwendeten Formeln in einem Dataframe gesammelt, was es auch möglich macht, die Anmerkungen "passend" untereinander zu positionieren. Vielen Leuten ist gar nicht bewußt, dass sich die einzelnen Layer eines Plots durchaus auf unterschiedliche Dataframes beziehen können!
Ich habe zwei in meinen Augen wesentliche Änderungen drin: zunächst habe ich den Plotdatensatz ins "lange Format" umgewandelt -
wodurch es möglich ist, verschiedene Linienzüge in unterschiedlichen Farben mit einem Statement abzufackeln.ggplot liebt das lange Datenformat (EDi - vor langer Zeit)
Zudem gewinnst Du die Freiheit, bei Bedarf die Grafiken im Handumdrehen über "facets" zu entflechten...
Zum zweiten habe ich die Anmerkungen zu den verwendeten Formeln in einem Dataframe gesammelt, was es auch möglich macht, die Anmerkungen "passend" untereinander zu positionieren. Vielen Leuten ist gar nicht bewußt, dass sich die einzelnen Layer eines Plots durchaus auf unterschiedliche Dataframes beziehen können!
Code: Alles auswählen
library(ggplot2)
library(hrbrthemes)
library(reshape)
#---Polynomiale Regression Einzelwerte aus Beispiel 1 Appendix C CSLI EP06-A----
#---Dateneingabe----
Sample <- data.frame(x = c(20, 20, 40, 40, 60, 60, 80, 80, 100, 100),
Messung = c(26.5, 26.2, 139, 138, 269, 273, 337, 343, 409, 404))
#---Berechnung Regressionsmodelle (linear, quadratisch "2" und kubisch "3")----
fit.1 <- lm(Messung ~ x, data=Sample)
Sample$linearPred <- predict(fit.1)
fit.2 <- lm(Messung ~ poly(x, 2, raw=TRUE), data=Sample)
Sample$quadraticPred <- predict(fit.2)
fit.3 <- lm(Messung ~ poly(x, 3, raw=TRUE), data=Sample)
Sample$cubicPred <- predict(fit.3)
Sample.long <- melt(Sample, id="x", variable_name="Reihe")
#---Formelerstellung lineare Regression für ggplot----
eqn.1 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.1)$coefficients[1,1], digits = 4),
b = format(summary(fit.1)$coefficients[2,1], digits = 4),
r2 = format(summary(fit.1)$r.squared, digits = 2)
))))
eqn.2 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)+c*italic(x)^2*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.2)$coefficients[1,1], digits = 4),
b = format(summary(fit.2)$coefficients[2,1], digits = 4),
c = format(summary(fit.2)$coefficients[3,1], digits = 4),
r2 = format(summary(fit.2)$r.squared, digits = 2)
))))
eqn.3 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)+c*italic(x)^2+d*italic(x)^3*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.3)$coefficients[1,1], digits = 4),
b = format(summary(fit.3)$coefficients[2,1], digits = 4),
c = format(summary(fit.3)$coefficients[3,1], digits = 4),
d = format(summary(fit.3)$coefficients[4,1], digits = 4),
r2 = format(summary(fit.3)$r.squared, digits = 2)
))))
Texte <- data.frame(lnr = 1:3, Reihe=c("linearPred","quadraticPred","cubicPred"), Equation = c(eqn.1, eqn.2, eqn.3))
#---Plot-Erstellung mit ggplot2---------
point <- ggplot(Sample.long) +
geom_point(aes(x, value, colour=Reihe)) +
geom_line(aes(x, value, colour=Reihe)) +
labs(title = 'Regressions Plot CLSI EP06-A Example 1 Appendix C',
subtitle = 'Polynomial regression Plot to test Linearity',
x = 'Concentration',
y = 'Ct-value',
caption = 'Andreas / andreas@random.com') +
geom_text(data=Texte, aes(x=25, y=440 - 20*lnr, label=(Equation), colour=Reihe), parse=TRUE, hjust="left") +
theme_ipsum() +
theme(plot.title = element_text(color = "#3a4fa7ff"),
plot.caption = element_text(color = "#3a4fa7ff",
face = 'bold')) +
scale_colour_manual(values=c("red","blue","darkgreen","blueviolet"))
point
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Polynomiale (polynomielle) Regression mit poly
Hi,
Ich werde jetzt mal weiter machen und das hoffentlich bald fertige dann mal zeigen.
Tja was soll ich sagen, genau dazu wollte ich was fragen, weil es in meinen plots nämlich besagte Verschiebungen gibt. Vielen Dank für dein Ausprobieren und damit für die Lösung zu meiner Angestellten Frage.
Ich werde jetzt mal weiter machen und das hoffentlich bald fertige dann mal zeigen.
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Polynomiale (polynomielle) Regression mit poly
danke für euere Antworten.
@Bigben:
bei dir wird ja eine Zahlenreihe zwischen 15 bis 405 in 5er Schritten erstellt und dann ein passendes Regressionsmodell gesucht ohne aber die Steigung der echten Regression zu berücksichtigen, oder stehe ich da auf dem Schlauch?!
Athomas macht es ja z.B. sobigben hat geschrieben: ↑Fr Apr 03, 2020 10:32 am Mein zweites Beispiel zeigt, wie man aus mit lm durchgeführte Regressionen Werte vorhersagen kann und wie man damit Kurven zeichnen kann. Das erfordert einen Umweg über lm und predict, dafür erhält man dann auch Werte, die man außerhalb von Abbildungen verwenden kann.
Code: Alles auswählen
fit.2 <- lm(Messung ~ poly(x, 2, raw=TRUE), data=Sample)
Sample$quadraticPred <- predict(fit.2)
Hier mal mein neuer Code angepasst (csv-Datei siehe Anhang) für ein echtes Beispiel. Ich muss log10 Skalierungen machen. Jetzte wäre es natürlich super wenn man die smooth-Kurven reinbekommt wie am Anfang evtl. mit Konfidenzintervallen. Es geht ja um den Nachweis von Linearität oder eben auch nicht.
Code: Alles auswählen
#---genutzte Bibliotheken----
library(ggplot2)
library(hrbrthemes)
library(reshape)
library(readr)
#---Einlesen Daten für Forum----
Sample1 <- read_csv("Sample1.csv", col_types = cols(x = col_number(),
y = col_number()))
#Notwendig weil sonst bei melt ein Fehler auftritt
x<- c(Sample1$x)
y<- c(Sample1$y)
Sample1<- data.frame(x, y)
#---Erster Plot zur Übersicht----
ggplot(Sample1, aes(x = x, y = y)) +
geom_point() + annotation_logticks(sides = "b") + scale_x_log10() +
geom_smooth(method="lm", color = 1) +
geom_smooth(method="lm", formula=y ~ I(x^2) + x, color = 2) +
geom_smooth(method="lm", formula=y ~ I(x^3) + I(x^2)+x, color = 3) +
theme_bw()
#---Berechnung Regressionsmodelle (linear, quadratisch "2" und kubisch "3")----
fit.1 <- lm(y ~ x, data=Sample1)
Sample1$linearPred <- predict(fit.1)
fit.2 <- lm(y ~ poly(x, 2, raw=TRUE), data=Sample1)
Sample1$quadraticPred <- predict(fit.2)
fit.3 <- lm(y ~ poly(x, 3, raw=TRUE), data=Sample1)
Sample1$cubicPred <- predict(fit.3)
Sample.long <- melt(Sample1, id="x", variable_name="Reihe")
#---Formelerstellung lineare Regression für ggplot----
eqn.1 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.1)$coefficients[1,1], digits = 4),
b = format(summary(fit.1)$coefficients[2,1], digits = 4),
r2 = format(summary(fit.1)$r.squared, digits = 4)
))))
eqn.2 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)+c*italic(x)^2*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.2)$coefficients[1,1], digits = 4),
b = format(summary(fit.2)$coefficients[2,1], digits = 4),
c = format(summary(fit.2)$coefficients[3,1], digits = 4),
r2 = format(summary(fit.2)$r.squared, digits = 4)
))))
eqn.3 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)+c*italic(x)^2+d*italic(x)^3*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.3)$coefficients[1,1], digits = 4),
b = format(summary(fit.3)$coefficients[2,1], digits = 4),
c = format(summary(fit.3)$coefficients[3,1], digits = 4),
d = format(summary(fit.3)$coefficients[4,1], digits = 4),
r2 = format(summary(fit.3)$r.squared, digits = 4)
))))
Texte <- data.frame(lnr = 1:3, Reihe=c("linearPred","quadraticPred","cubicPred"), Equation = c(eqn.1, eqn.2, eqn.3))
point <- ggplot(Sample.long) +
annotation_logticks(sides = "b") + scale_x_log10() +
geom_point(aes(x, value, colour=Reihe)) +
geom_line(aes(x, value, colour=Reihe)) +
labs(title = 'Regressions Plot CLSI EP06-A Example 1 Appendix C',
subtitle = 'Polynomial regression Plot to test Linearity',
x = 'Concentration',
y = 'Ct-value',
caption = 'Andreas / andreas@random.com') +
geom_text(data=Texte, aes(x=25, y=40 - 1*lnr, label=(Equation), colour=Reihe), parse=TRUE, hjust="left") +
theme_ipsum() +
theme(plot.title = element_text(color = "#3a4fa7ff"),
plot.caption = element_text(color = "#3a4fa7ff",
face = 'bold')) +
scale_colour_manual(values=c("red","blue","darkgreen","blueviolet"))
point
Ich hoffe es ist nicht ganz unverständlich was ich geschrieben habe.
Disclaimer :
Nur zur Info ich bin kein Statistiker sondern Biologe, der aber alles können muss laut Chef!!! Ihr dürft gerne meine unstatistischen Aussagen korrigieren. Ich wäre sogar dankbar.
Ihr macht hier eine richtig tolle Arbeit.
Ich mag R und dank euch lerne ich langsam damit umzugehen!
Zuletzt geändert von Werekorden am Di Apr 07, 2020 11:58 am, insgesamt 1-mal geändert.
Re: Polynomiale (polynomielle) Regression mit poly
Ich guck mir das Zeug mal an, aber zunächst: diese Irrsinns-Liste von nicht benötigten Packages geht mir tierisch auf den Senkel - das irritiert mehr, als dass es hilft !
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Polynomiale (polynomielle) Regression mit poly
Hi,
Ich bin ein Depp. Ich habe vergessen, dass die Lm-Methode ja die Info braucht, dass es sich um Log Werte handelt. Sorry habe ich in meinem Code vergessen bzw. noch nicht ausprobiert gehabt.
Ok die Liste habe ich drin, weil ich die später mit Cairo usw. exportieren wollte und gestehe ich habe vergessen welche notwendig waren.
Sind raus:
Ich bin ein Depp. Ich habe vergessen, dass die Lm-Methode ja die Info braucht, dass es sich um Log Werte handelt. Sorry habe ich in meinem Code vergessen bzw. noch nicht ausprobiert gehabt.
Ok die Liste habe ich drin, weil ich die später mit Cairo usw. exportieren wollte und gestehe ich habe vergessen welche notwendig waren.
Sind raus:
Zuletzt geändert von Werekorden am Di Apr 07, 2020 5:40 pm, insgesamt 1-mal geändert.
Re: Polynomiale (polynomielle) Regression mit poly
Ich verkneife mir mal diverse Kommentare zu Deinem Programm, weil ich befürchte, dass die eher verwirren würden.Ich habe vergessen, dass die Lm-Methode ja die Info braucht, dass es sich um Log10 Werte handelt.
Haben wir jetzt fertig oder bleibt da noch was ?
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Polynomiale (polynomielle) Regression mit poly
Hallo Athomas,
Gerne Kommentare. Ich will ja lernen.
Ich denke das war es . Ich hatte halt zuerst das Beispiel der CLSI genutzt um es dann für meinen spezifische Aufgabe anzupassen.
Nochmals danke bis hierher für eure tolle Hilfe. Ich merke ich bin gerade ein wenig auf.
Hier noch ein alter Code wo ich mal nur Lineare Regression für Log Daten gemacht hatte:
VG
Andreas
Gerne Kommentare. Ich will ja lernen.
Ich denke das war es . Ich hatte halt zuerst das Beispiel der CLSI genutzt um es dann für meinen spezifische Aufgabe anzupassen.
Nochmals danke bis hierher für eure tolle Hilfe. Ich merke ich bin gerade ein wenig auf.
Hier noch ein alter Code wo ich mal nur Lineare Regression für Log Daten gemacht hatte:
Code: Alles auswählen
#------CMV-Linearität mit ggplot erstellt----------------------------------------
library(ggplot2)
library(hrbrthemes)
#Dateneingabe
Konz_x <- c(10000, 10000, 10000, 1000, 1000, 1000, 100, 100, 100)
CT_Wert_y <- c(21.9053, 21.7273, 21.6431, 25.4452, 25.1422, 25.4432, 28.5071, 28.6153, 28.4508)
data.frame(Konz_x, CT_Wert_y)
p <- data.frame(Konz_x, CT_Wert_y)
#Model-Erstellung für Formeldarstellung in plot
model <- lm(CT_Wert_y ~ log10(Konz_x))
summary(model)
model
#Formel-Erstellung
eqn <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(model)$coefficients[1,1], digits = 4),
b = format(summary(model)$coefficients[2,1], digits = 4),
r2 = format(summary(model)$r.squared, digits = 5)
))))
eqn
parse(text = eqn)
#Plot-Erstellung
point <- ggplot(p, aes(x=Konz_x, y=CT_Wert_y))
p <- point + geom_point() + annotation_logticks(sides = "b") + scale_x_log10() + stat_smooth(method=lm, se=0.95, fullrange = TRUE, colour="purple") +
labs(title = 'Linearity Plot of CMV-Kit',
subtitle = 'Linearity Plot to find Range of Linearity with Coefficient of Determination',
x = 'Copies/ml patient sample',
y = 'Ct-value',
caption = 'Dr. Andreas Breß / bress@andiatec.com') +
annotate("text", label = eqn, parse = TRUE, x=2600 , y=20.5) +
theme_ipsum() +
theme(plot.title = element_text(color = "#3a4fa7ff"), plot.caption = element_text(color = "#3a4fa7ff", face = 'bold'))
p
Andreas
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Polynomiale (polynomielle) Regression mit poly
Hallo Athomas und Bigben
Ich komme irgendwie nicht weiter. Der Y-intercepter müsste in diesem Beispiel irgendwo bei y ≈ 40 liegen. Laut linearer Regressions-Funktion liegt er aber bei ≈ 23. Die Steigung passt mit 3.32
Irgendwie muss es damit zusammenhängen das ein logarithmus von kleiner 1 eine negative Zahl ausgibt, aber irgendwie ist bei mir ein Knoten im Hirn.
Mit den auskommentierten Werten klappt es aber eben nicht mit Werten unter 1.
Es wäre wirklich toll wenn Ihr mir helfen könntet.
VG
Andreas
Ich komme irgendwie nicht weiter. Der Y-intercepter müsste in diesem Beispiel irgendwo bei y ≈ 40 liegen. Laut linearer Regressions-Funktion liegt er aber bei ≈ 23. Die Steigung passt mit 3.32
Irgendwie muss es damit zusammenhängen das ein logarithmus von kleiner 1 eine negative Zahl ausgibt, aber irgendwie ist bei mir ein Knoten im Hirn.
Mit den auskommentierten Werten klappt es aber eben nicht mit Werten unter 1.
Code: Alles auswählen
#---genutzte Libraries----
library(ggplot2)
library(hrbrthemes)
library(reshape)
#---Dateneingabe für jede 10er Faktor in x wird y um 3.32 größer----
x<-c(0.1, 0.1, 0.1, 0.01, 0.01, 0.01, 0.001, 0.001, 0.001, 0.0001, 0.0001, 0.0001)
y<-c(26.32, 26.32, 26.32, 29.64, 29.64, 29.64, 32.96, 32.96, 32.96, 36.28, 36.28, 36.28)
#x<-c(10000, 10000, 10000, 1000, 1000, 1000, 100, 100, 100, 10, 10, 10)
#y<-c(26.32, 26.32, 26.32, 29.64, 29.64, 29.64, 32.96, 32.96, 32.96, 36.28, 36.28, 36.28)
sample1 <- data.frame(x, y)
#---Berechnung Regressionsmodelle (linear, quadratisch "2" und kubisch "3") für Formel-Erstellung----
fit.1 <- lm(y ~ (log10(x)), data=sample1)
fit.2 <- lm(y ~ poly(log10(x), 2, raw=TRUE), data=sample1)
fit.3 <- lm(y ~ poly(log10(x), 3, raw=TRUE), data=sample1)
#---Formelerstellung lineare Regression für ggplot----
eqn.1 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.1)$coefficients[1,1], digits = 4),
b = format(summary(fit.1)$coefficients[2,1], digits = 4),
r2 = format(summary(fit.1)$r.squared, digits = 4)
))))
eqn.2 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)+c*italic(x)^2*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.2)$coefficients[1,1], digits = 4),
b = format(summary(fit.2)$coefficients[2,1], digits = 4),
c = format(summary(fit.2)$coefficients[3,1], digits = 4),
r2 = format(summary(fit.2)$r.squared, digits = 4)
))))
eqn.3 <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)+c*italic(x)^2+d*italic(x)^3*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit.3)$coefficients[1,1], digits = 4),
b = format(summary(fit.3)$coefficients[2,1], digits = 4),
c = format(summary(fit.3)$coefficients[3,1], digits = 4),
d = format(summary(fit.3)$coefficients[4,1], digits = 4),
r2 = format(summary(fit.3)$r.squared, digits = 4)
))))
Texte <- data.frame(lnr = 1:3, Reihe=c("linearPred","quadraticPred","cubicPred"), Equation = c(eqn.1, eqn.2, eqn.3))
#---Plot-Erstellung----
plot<- ggplot(sample1, aes(x = x, y = y)) +
annotation_logticks(sides = "b") + scale_x_log10() +
labs(title = 'Regressions Plot PCR Results',
subtitle = 'Polynomial regression Plot to test Linearity',
x = 'Concentration IU/ml',
y = 'Cq-value',
caption = 'Andreas / andreas@random.com')+
geom_point(aes(x=x, y=y)) +
geom_smooth(method="lm", color = 3, fullrange = TRUE) +
geom_smooth(method="lm", formula=y ~ I(x^2) + x, color = 4, fullrange = TRUE) +
geom_smooth(method="lm", formula=y ~ I(x^3) + I(x^2) + x, color = 2, fullrange = TRUE) +
geom_text(data=Texte, aes(x=0.0001, y=28 - 0.75*lnr, label=(Equation), colour=Reihe), parse=TRUE, hjust="left",
show.legend=FALSE) +
theme_bw()
plot
VG
Andreas
Re: Polynomiale (polynomielle) Regression mit poly
Also wenn ich Deinen Code mal auf das nötige zusammenstreiche und Deine Punkte plotte, dann passt der Intercept von 23 ziemlich gut
LG,
Bernhard
Code: Alles auswählen
library(ggplot2)
x <- c(0.1, 0.1, 0.1, 0.01, 0.01, 0.01, 0.001, 0.001, 0.001, 0.0001, 0.0001, 0.0001)
y <- c(26.32, 26.32, 26.32, 29.64, 29.64, 29.64, 32.96, 32.96, 32.96, 36.28, 36.28, 36.28)
sample1 <- data.frame(x, y)
fit.1 <- lm(y ~ log10(x), data=sample1)
summary(fit.1)
ggplot(sample1, aes(x=log10(x), y=y)) +
geom_count() +
geom_abline(slope= -3.32, intercept = 23) +
geom_point(aes(x=0, y=23),pch = "X", size = 20, color="orange")
Bernhard
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte