Konfidenzintervall eines Punktes in Regressionsgeraden

Modelle zur Korrelations- und Regressionsanalyse

Moderator: EDi

Antworten
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Konfidenzintervall eines Punktes in Regressionsgeraden

Beitrag von Werekorden »

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
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Konfidenzintervall eines Punktes in Regressionsgeraden

Beitrag von EDi »

Code: Alles auswählen

if (F) {  
:shock:
Kann man machen - sollte man aber nicht. if (F) { ist schon hart zu sehen... (unit tests sind ja sinnvoll, aber so???)
Hat da jemand eine Idee wie ich das aus "upr" und "lwr" rausziehen kann.
Die funktion gibt das zurück,

Code: Alles auswählen

 invisible(list(r=r, lod95=lod95))
invisble heißt nur, das die print-methode nicht aufgerufen wird. Es wird also eine Liste mit r und lod95 zurückgegeben.
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
Bild.
Werekorden
Beiträge: 83
Registriert: So Feb 04, 2018 7:52 pm

Re: Konfidenzintervall eines Punktes in Regressionsgeraden

Beitrag von Werekorden »

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

Re: Konfidenzintervall eines Punktes in Regressionsgeraden

Beitrag von Werekorden »

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?
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Konfidenzintervall eines Punktes in Regressionsgeraden

Beitrag von EDi »

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:

Code: Alles auswählen

grid <- seq(0.00001, max(1, dat$conc), length.out=10000)
Besser wäre es, ein log-grid zu nehmen, anstatt das gleichmäßig hier.

2. Schau dir mal diesen plot an:

Code: Alles auswählen

plot(grid, fit, log = 'x')
lines(grid, upr)
lines(grid, lwr)
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):

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

Re: Konfidenzintervall eines Punktes in Regressionsgeraden

Beitrag von Werekorden »

Hallo Edi,


wow danke für die Hilfe soweit. :P

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))])
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
Antworten