Konfidenzintervall eines Punktes in Regressionsgeraden
Moderator: EDi
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Konfidenzintervall eines Punktes in Regressionsgeraden
Hi,
Ich sollte folgendes R-Skript nutzen und bin nicht so firm mit R.
http://frisby5.blogspot.com/2016/02/pro ... ction.html
Jetzt benötige ich die Werte für das Konfidenzintervall bei meinem LOD-Punkt.
Hat da jemand eine Idee wie ich das aus "upr" und "lwr" rausziehen kann.
Ich könnte auch diese Skript nutzen wobei für mich das erste wesentlich einfacher zu verstehen ist:
https ://github.com/cmerkes/qPCR_LOD_Calc, and our additional analysis code can be found at https ://github.com/cmerkes/LOD_Analysis. The data used in this study are available at https ://doi.org/10.5066/P9AKHU1R
Das ganze stammt aus diesem, Paper:
Klymus KE, Merkes CM, Allison MJ,et al. Reporting the limits of detection and quantification forenvironmental DNA assays. Environmental DNA. 2019;00: 1–1 2 . https ://doi.org/10.1002/edn3.29
Dank euch
Andreas
Ich sollte folgendes R-Skript nutzen und bin nicht so firm mit R.
http://frisby5.blogspot.com/2016/02/pro ... ction.html
Jetzt benötige ich die Werte für das Konfidenzintervall bei meinem LOD-Punkt.
Hat da jemand eine Idee wie ich das aus "upr" und "lwr" rausziehen kann.
Ich könnte auch diese Skript nutzen wobei für mich das erste wesentlich einfacher zu verstehen ist:
https ://github.com/cmerkes/qPCR_LOD_Calc, and our additional analysis code can be found at https ://github.com/cmerkes/LOD_Analysis. The data used in this study are available at https ://doi.org/10.5066/P9AKHU1R
Das ganze stammt aus diesem, Paper:
Klymus KE, Merkes CM, Allison MJ,et al. Reporting the limits of detection and quantification forenvironmental DNA assays. Environmental DNA. 2019;00: 1–1 2 . https ://doi.org/10.1002/edn3.29
Dank euch
Andreas
Re: Konfidenzintervall eines Punktes in Regressionsgeraden
Code: Alles auswählen
if (F) {
Kann man machen - sollte man aber nicht. if (F) { ist schon hart zu sehen... (unit tests sind ja sinnvoll, aber so???)
Die funktion gibt das zurück,Hat da jemand eine Idee wie ich das aus "upr" und "lwr" rausziehen kann.
Code: Alles auswählen
invisible(list(r=r, lod95=lod95))
Pass doch den code an damit dir upr & lwr auch zurückgegeben wird...
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.
Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
.
Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
.
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Konfidenzintervall eines Punktes in Regressionsgeraden
Hi Edi,
Die Funktion von der Homepage ist wie gesagt nicht von mir. Ich versuche das nur zu verstehen und bin sehr DAUig drauf was solche komplexeren Funktionen angeht. Werde später mal schauen ob ich mit deinen tIps weiterkomme.
Danke soweit.
Andreas
Die Funktion von der Homepage ist wie gesagt nicht von mir. Ich versuche das nur zu verstehen und bin sehr DAUig drauf was solche komplexeren Funktionen angeht. Werde später mal schauen ob ich mit deinen tIps weiterkomme.
Danke soweit.
Andreas
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Konfidenzintervall eines Punktes in Regressionsgeraden
Hi Edi,
Ok das mir upr und lwr die Zahlenwerte rausschmeisst ist mir schon klar, nur verstehe ich diese nicht. Ich kann nicht erkennen was das für Werte sind?
Deshalb weiß ich auch nicht wie ich den passenden Wert für meinen LoD-Wert finden kann.
Ich stehe gerade leider total auf dem Schlauch.
ODER:
bedeuten die Werte in "lwr, upr" und "fit" etwa die Prozentangaben für bestimmte Werte der x-Achse, also für die genome/ml sample und das liegt daran das für "lwr, upr, fit" die das inverse der Link-Funktion genutzt wurde?
Ok das mir upr und lwr die Zahlenwerte rausschmeisst ist mir schon klar, nur verstehe ich diese nicht. Ich kann nicht erkennen was das für Werte sind?
Deshalb weiß ich auch nicht wie ich den passenden Wert für meinen LoD-Wert finden kann.
Ich stehe gerade leider total auf dem Schlauch.
ODER:
bedeuten die Werte in "lwr, upr" und "fit" etwa die Prozentangaben für bestimmte Werte der x-Achse, also für die genome/ml sample und das liegt daran das für "lwr, upr, fit" die das inverse der Link-Funktion genutzt wurde?
Re: Konfidenzintervall eines Punktes in Regressionsgeraden
Du hast den y- Wert bei 0.95.
Die Konfidenzbänder sind in upr & lwr gegeben.
Du kannst jetzt jetzt einfach nur schauen wo upw & lwr die 0.95 kreuzt . Es gibt sichlich auch geschlossene Lösung, aber die kannich jetzt heute Abend nicht herleiten
1. verfeinere dein Grid:
Besser wäre es, ein log-grid zu nehmen, anstatt das gleichmäßig hier.
2. Schau dir mal diesen plot an:
Das erklärt was wir haben wollen. Die Schnittpunkte der 3 Linie mit der horizontalen bei 0.95
3. Einfach nur schauen wo die 0.95 gekreuzt wird (dazu muss dein grid fein genug sein):
Ich würde dafür drc verwenden...
https://www.crcpress.com/Dose-Response- ... 1138034310
Hier ein Beispiel (muss nicht zu deinem Problem passen):
Die Konfidenzbänder sind in upr & lwr gegeben.
Du kannst jetzt jetzt einfach nur schauen wo upw & lwr die 0.95 kreuzt . Es gibt sichlich auch geschlossene Lösung, aber die kannich jetzt heute Abend nicht herleiten
1. verfeinere dein Grid:
Code: Alles auswählen
grid <- seq(0.00001, max(1, dat$conc), length.out=10000)
2. Schau dir mal diesen plot an:
Code: Alles auswählen
plot(grid, fit, log = 'x')
lines(grid, upr)
lines(grid, lwr)
3. Einfach nur schauen wo die 0.95 gekreuzt wird (dazu muss dein grid fein genug sein):
Code: Alles auswählen
grid[max(which(lwr >= 0.95))]
grid[max(which(upr > 0.95))]
Ich würde dafür drc verwenden...
https://www.crcpress.com/Dose-Response- ... 1138034310
Hier ein Beispiel (muss nicht zu deinem Problem passen):
Code: Alles auswählen
library(drc)
head(earthworms)
mod <- drm(number/total ~ dose,
weights = total,
data = earthworms,
fct = LN.2(),
type = "binomial")
plot(mod, type = 'confidence', yt = seq(0, 1, 0.1))
plot(mod, add = TRUE, type = 'obs')
# EC10, 50 & 90 + CI
eds <- ED(mod, c(10, 50, 90), interval = 'delta')
abline(v = eds[2, 1])
abline(h = 0.5, lty = 'dotted')
Bitte immer ein reproduzierbares Minimalbeispiel angeben. Meinungen gehören mir und geben nicht die meines Brötchengebers wieder.
Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
.
Dieser Beitrag ist lizensiert unter einer CC BY 4.0 Lizenz
.
-
- Beiträge: 83
- Registriert: So Feb 04, 2018 7:52 pm
Re: Konfidenzintervall eines Punktes in Regressionsgeraden
Hallo Edi,
wow danke für die Hilfe soweit.
Ich hab das mal in "meinen" Code eingefügt:
Seltsamerweise, zeigt er aber den Plot-Titel und die Achsenbeschriftung nur an wenn man die Schritte einzeln ausführt und nicht wenn man die Funktion durchlaufen lässt.
Außerdem gibt dein:
jeweils den Wert 1 aus, was ja nicht so ganz stimmen kann oder habe ich da wieder einen Denkfehler?
Auch hier muss man es übrigens einzeln ausführen und es funktioniert nicht beim ausführen der Funktion.
Trotzdem schon einmal vielen Dank bis hierher.
Andreas
wow danke für die Hilfe soweit.
Ich hab das mal in "meinen" Code eingefügt:
Code: Alles auswählen
lod.model <- function(dat)
{
## PURPOSE: Make a Probit Model Fitting and Ploting for LOD Data
##
## ARGUMENT:
## dat: a data frame with three columns "conc", "success", "failure"
## "conc" is the Target Concentration Levels.
## "success" is the Number of Positive Results.
## "failure" is the Number of Negative Results.
##
## RETURN: Plot, GLM model object, and 95% LOD Estimate.
##
## AUTHOR: Feiming Chen, Date: 10 Feb 2016, 14:28
## Per "CLSI EP17A2E" Guidance, apply default recommendation:
## 1. Remove Zero Predictor Values.
## 2. Appply "log10" transformation on predictors.
## 3. Use Probit Model
## 4. Report 95% LOD (concentration level)
## remove 0 levels in "conc" since will do "log10" transform on predictor.
dat <- subset(dat,)
r <- glm(cbind(success, failure) ~ log10(conc), data=dat, family=binomial(link="probit"))
grid <- seq(0.00001, max(1, dat$conc), length.out=10000)
preds <- predict(r, newdata=data.frame(conc=grid), type="link", se.fit=TRUE)
se <- 1.96 * preds$se.fit
upr <- unlist(r$family$linkinv(preds$fit + se))
lwr <- unlist(r$family$linkinv(preds$fit - se))
fit <- unlist(r$family$linkinv(preds$fit))
ce <- coef(r)
lod95 <- 10^((r$family$linkfun(0.95) - ce[1])/ce[2]) # 95% LOD Estimate
## Plot pch=19 bedeutet solid circles werden gemalt, log="X" gibt an, dass die X-achse im logscale sein soll
# with(dat, plot(conc, success/(success+failure),ylab="Probability", xlab="Concentration", log="x", pch=1, ylim=c(0, 1)))
#title(main=paste("Probit Model, 95% LOD =", round(lod95, 3)))
# lines(grid, fit)
# lines(grid, lwr, lty=2)
# lines(grid, upr, lty=2)
# abline(h=0.95, lty=2, col="red")
# abline(v=lod95, lty=2, col="red")
# axis(2, 0.95, labels = "0.95")
plot(grid, fit, log = 'x', ylab="Probability", xlab="Concentration", pch=1)
title(main=paste("Probit Model, 95% LOD =", round(lod95, 3)))
lines(grid, upr, lty=2)
lines(grid, lwr, lty=2)
abline(h=0.95, lty=2, col="red")
abline(v=lod95, lty=2, col="red")
axis(2, 0.95, labels = "0.95")
(grid[max(which(lwr >= 0.95))])
(grid[max(which(upr > 0.95))])
invisible(list(r=r, lod95=lod95))
}
if (F) { # Unit Test
dat <- data.frame(conc=c(5.4e5, 1.1e5, 5e4, 1e4, 5e3, 1e3),
success=c(6, 6, 6, 3, 0, 0),
failure=c(0, 0, 0, 3, 6, 6))
lod.model(dat)
}
## CLSI Example: EP17A2E Table C1 & C2, Lot 1, LOD=0.077 CFU/mL.
dat <- data.frame(conc=c(0.025, 0.05, 0.15, 0.3, 0.5, 0.006, 0.014),
success=c(23, 29, 32, 32, 32, 11, 15),
failure=c(9, 3, 0, 0, 0, 19, 15))
r <- lod.model(dat) # Verify 95% LOD = 0.077
Seltsamerweise, zeigt er aber den Plot-Titel und die Achsenbeschriftung nur an wenn man die Schritte einzeln ausführt und nicht wenn man die Funktion durchlaufen lässt.
Außerdem gibt dein:
Code: Alles auswählen
(grid[max(which(lwr >= 0.95))])
(grid[max(which(upr > 0.95))])
Auch hier muss man es übrigens einzeln ausführen und es funktioniert nicht beim ausführen der Funktion.
Trotzdem schon einmal vielen Dank bis hierher.
Andreas