Polynomiale (polynomielle) Regression mit poly

Modelle zur Korrelations- und Regressionsanalyse

Moderator: EDi

bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von bigben »

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
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Athomas »

So, jetzt habe ich mich auch noch mal am Aufhübschen Deiner Plots versucht - Hausarrest macht's möglich :D !
Ich habe zwei in meinen Augen wesentliche Änderungen drin: zunächst habe ich den Plotdatensatz ins "lange Format" umgewandelt -
ggplot liebt das lange Datenformat (EDi - vor langer Zeit)
wodurch es möglich ist, verschiedene Linienzüge in unterschiedlichen Farben mit einem Statement abzufackeln.
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
Dateianhänge
Linienzüge.jpeg
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Werekorden »

Hi,
Athomas hat geschrieben: Fr Apr 03, 2020 11:31 amVielen Leuten ist gar nicht bewußt, dass sich die einzelnen Layer eines Plots durchaus auf unterschiedliche Dataframes beziehen können!
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.
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Werekorden »

Sample1.csv
Sample1
(834 Bytes) 46-mal heruntergeladen
Hallo Bigben, Hallo Athomas,

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?!
bigben 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.
Athomas macht es ja z.B. so

Code: Alles auswählen

fit.2                <- lm(Messung ~ poly(x, 2, raw=TRUE), data=Sample)
Sample$quadraticPred <- predict(fit.2)
Damit wird zunächst angegeben welche Formel-Typ man nimmt und mit welchen Daten dann das Regressionsmodel errechnet wird. Danach dann per predict werden Werte basierend auf dem Model erstellt und in die Tabelle eingefügt.



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.
Sample1.csv
Sample1
(834 Bytes) 46-mal heruntergeladen

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. :P
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.
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Athomas »

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 :evil: !
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Werekorden »

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: :|
Zuletzt geändert von Werekorden am Di Apr 07, 2020 5:40 pm, insgesamt 1-mal geändert.
Athomas
Beiträge: 768
Registriert: Mo Feb 26, 2018 8:19 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Athomas »

Ich habe vergessen, dass die Lm-Methode ja die Info braucht, dass es sich um Log10 Werte handelt.
Ich verkneife mir mal diverse Kommentare zu Deinem Programm, weil ich befürchte, dass die eher verwirren würden.
Haben wir jetzt fertig oder bleibt da noch was :D ?
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Werekorden »

Hallo Athomas,

Gerne Kommentare. Ich will ja lernen.
Ich denke das war es :oops: . 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
VG
Andreas
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von Werekorden »

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.

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
Es wäre wirklich toll wenn Ihr mir helfen könntet.

VG
Andreas
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: Polynomiale (polynomielle) Regression mit poly

Beitrag von bigben »

Also wenn ich Deinen Code mal auf das nötige zusammenstreiche und Deine Punkte plotte, dann passt der Intercept von 23 ziemlich gut

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")
LG,
Bernhard
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Antworten