Hi,
ja das habe ich mittlerweile auch herausgefunden. Ich versuche gerade ein CLSI-Beispiel nachzustellen damit wir das dann nutzen können für unserer Linearitäts-Austestungen.
Leider scheint es einen Fehler in dem CLSi-Dokument zu geben. Die bekommen ein anderes b1 raus als R bei der BErechnung der lm ersten Grades.
Ich habs mal in ggplot gemacht mit Formel
aber jetzt steh ich wieder auf dem Schlauch mit dem Überlagern der Regressionsmodellen zweiten und dritten Grades.
Code: Alles auswählen
library(Rcpp)
library(rlang)
library(gcookbook)
library(ggplot2)
library(hrbrthemes)
library(magick)
library(plyr)
library(dplyr)
library(MASS)
library(scales)
library(Cairo)
library(extrafont)
x<-c(20, 20, 40, 40, 60, 60, 80, 80, 100, 100)
y<-c(26.5, 26.2, 139, 138, 269, 273, 337, 343, 409, 404)
sample1 <- data.frame(y, x)
sample1
plot(sample1$x, sample1$y, type = "l")
fit1<- lm(sample1$y ~ sample1$x)
summary(fit1)
fit2<- lm(sample1$y ~ poly(sample1$x, 2, raw=TRUE))
summary(fit2)
fit3<- lm(sample1$y ~ poly(sample1$x, 3, raw=TRUE))
summary(fit3)
plot(sample1$x, sample1$y, type = "l", lwd=1, main = "Vergleich Polynomiale Regression", xlab = "Konzentration", ylab = "Wert")
points(sample1$x, predict(fit1), type="l", lty=2, col="black", lwd=1)
points(sample1$x, predict(fit2), type="l", lty=1, col="red", lwd=2)
points(sample1$x, predict(fit3), type="l", lty=3, col="blue", lwd=2)
#--------Formelerstellung------------
eqn <- as.character(as.expression(substitute(italic(y)==a+b*italic(x)*"," ~~ italic(r)^2~ "=" ~ r2,
list(a = format(summary(fit1)$coefficients[1,1], digits = 4),
b = format(summary(fit1)$coefficients[2,1], digits = 4),
r2 = format(summary(fit1)$r.squared, digits = 2)
))))
eqn
parse(text = eqn)
#---------ggPlot-Erstellung----------------
point <- ggplot(sample1, aes(x=x, y=y))
point + geom_point() + stat_smooth(method=lm, se=0.99, fullrange = TRUE, colour="purple") +
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') +
annotate("text", label = eqn, parse = TRUE, x=50 , y=400) +
theme_ipsum() +
theme(plot.title = element_text(color = "#3a4fa7ff"), plot.caption = element_text(color = "#3a4fa7ff", face = 'bold'))