Beta-Funktion auf Methylomdaten fitten

Modelle zur Korrelations- und Regressionsanalyse

Moderator: EDi

Antworten
Curnen
Beiträge: 27
Registriert: Fr Nov 18, 2016 3:45 pm

Beta-Funktion auf Methylomdaten fitten

Beitrag von Curnen »

Hallo liebe Forum-Helfer,

Ich bräuchte Hilfe dabei, die Beta-Funktion auf meine Daten zu fitten. Es handelt sich dabei um Messdaten aus epigenetischen Analysen, die den Wertebereich zwischen 0 und 1 annehmen können. Ihre Häufigkeitsverteilung wird in etlichen Fachpublikationen mit der Beta-Funktion beschrieben (wobei shape1=shape2 üblicherweise vereinfachend angenommen wird). Für einzelne Abschnitte des Genoms gibt es Referenzwerte für den Alpha-Wert=shape1=shape2 und ich würde für meine Daten gerne schauen, wie nah ich an diesen Referenzwerten liege.

Zur Veranschaulichung hier als Plot:

Code: Alles auswählen

methscores <-  sample(seq(0,1,by=0.0625),1e4, replace = TRUE, prob = c(40,1,2,1,4,1,2,1,8,1,2,1,4,1,2,1,100))

# calculate density
methscores_dens <- density(methscores)
plotdata <- data.frame("x"=methscores_dens[["x"]],"y"=methscores_dens[["y"]])

library("ggplot2")

ggplot(plotdata, aes(x, y)) + geom_line() + 
  geom_line(data=subset(plotdata,x>=0 & x<=1),color="blue",size=1.1) +
  stat_smooth(method = "lm", formula = y ~ dbeta(x,shape1=0.5,shape2=0.5),color="red",se=FALSE)


Der blau hervorgehobene Bereich von 0 bis 1 soll durch eine Beta-Funktion (in rot) möglichst gut beschrieben werden. Leider sehen die Messdaten tatsächlich so asymmetrisch aus wie im Beispiel - egal ob ich die echten Daten oder die simulierten aus dem Beispiel nehme, erhalte ich diese Fehlermeldung:

Code: Alles auswählen

library("MASS")
fit.beta <- fitdistr(methscores,"beta", start = list(shape1=0.5, shape2=0.5))
Error in stats::optim(x = c(0, 1, 0.1875, 0.375, 1, 1, 1, 1, 1, 1, 1,  : 
  initial value in 'vmmin' is not finite


Mit wirklich aus einer Beta-Funktion stammenden Werten funktioniert es so (zumindest wenn man die Warnungen außen vor lässt :? )

Code: Alles auswählen

library("MASS")
fitdistr(rbeta(1e4,shape1=0.7, shape2=0.4),"beta", start = list(shape1=0.5, shape2=0.5))
  shape1        shape2   
  0.710951827   0.400015091 
 (0.009944133) (0.004775643)
Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs wurden erzeugt
2: In densfun(x, parm[1], parm[2], ...) : NaNs wurden erzeugt
[...]


Habt ihr eine Idee, wie ich trotzdem eine Beta-Funktion gefittet bekomme, selbst wenn sie die zugrundeliegenden Daten eher schlecht beschreiben dürfte?
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von bigben »

Hi!
100 indicates that optim encountered an internal error.
Nicht sehr hilfreich.

Ich kenn mich da nicht aus, aber wenn es Dir nicht so ganz genau auf die Methode ankommt, habe ich hier wenigstens einen Weg zu einem Ergebnis:

Code: Alles auswählen

> library(fitdistrplus)
> fitdist(methscores, "beta", start=list(shape1=.1, shape2=.1))
<simpleError in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data,     gr = gradient, ddistnam = ddistname, hessian = TRUE, method = meth,     lower = lower, upper = upper, ...): function cannot be evaluated at initial parameters>
Error in fitdist(methscores, "beta", start = list(shape1 = 0.1, shape2 = 0.1)) : 
  the function mle failed to estimate the parameters, 
                with the error code 100
> fitdist(methscores, "beta", method="mme", start=list(shape1=.1, shape2=.1))
Fitting of the distribution ' beta ' by matching moments 
Parameters:
         estimate
shape1 0.13602291
shape2 0.06636766
LG,
Bernhard


EDIT: Sorry, aber das mit dem mme ist vielleicht nicht so schlau. Visuell geschätzt müsste die Kurve rechts bei der 1 höher sein:

Code: Alles auswählen

> hist(methscores, probability = TRUE)
> curve(dbeta(x, 0.136, 0.606), add = TRUE)
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von EDi »

Die beta-Verteilung schließt die 0 und 1 nicht mit ein - könnte das Problem erklären (mir scheint du hast viele solcher Werte).
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.
Curnen
Beiträge: 27
Registriert: Fr Nov 18, 2016 3:45 pm

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von Curnen »

Hallo ihr Beiden,

Danke für eure Zeit und die Antworten.
EDi hat geschrieben: Di Nov 14, 2017 7:32 pm Die beta-Verteilung schließt die 0 und 1 nicht mit ein - könnte das Problem erklären (mir scheint du hast viele solcher Werte).


Jep, das war tatsächlich ein Problem, das ich mir ungefähr zeitgleich mit deiner Antwort dann endlich auch zusammengegoogelt hatte.

Code: Alles auswählen

methscores <-  sample(seq(0,1,by=0.0625),1e4, replace = TRUE, prob = c(40,1,2,1,4,1,2,1,8,1,2,1,4,1,2,1,100))
methylscores[methylscores<=0] <- 0.0001
methylscores[methylscores>=1] <- 0.9999
So richtig zuverlässig läuft das leider immer noch nicht - für manche Subsamples der Messdaten funktioniert es schon für andere noch nicht:

Code: Alles auswählen

Error in stats::optim(x = c(0.9091, 0.8889, 0.9, 0.9231, 0.9091, 0.4545,  : 
  non-finite finite-difference value [2]
Ich werde daher jetzt mal mit dem Paket von bigben experimentieren.Würdet ihr denn empfehlen rbeta auf die originalen Werte oder pbeta auf die y-Werte der Density zu fitten?

Fun fact am Rande: Die Kumaraswamy-Verteilung - von der ich bislang noch nie was gehört hatte, aber die wohl eine leichter zu handhabende Alternative zur Beta-Verteilung sei - fittet R schon jetzt anstandslos:
kumaraswamy <- function(x,a,b){a*b*x^(a-1)*(1-x^a)^(b-1)}
bigben
Beiträge: 2771
Registriert: Mi Okt 12, 2016 9:09 am

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von bigben »

Woran der EDi wieder alles denkt. Wikipedia schreibt
Die Betaverteilung ist eine stetige Wahrscheinlichkeitsverteilung über dem Intervall (0,1). [...] Für p,q >=1 lässt sich (0,1) durch [0,1] ersetzen.
Leider haben wir es hier mit p und q kleiner 1 zu tun. Ich dachte eben, man könne die Daten vielleicht ein wenig stauchen, dass Sie in (0,1) hinein passen:

Code: Alles auswählen

meth <- methscores * .99998 + .00001
Dann gelingt zwar das Fitten, aber die Parameter verändern sich schon recht deutlich, wenn man mehr oder weniger Nachkommastellen einsetzt. Könnte man trotzdem weiter verfolgen, falls es bei sehr vielen Nachkommastellen wieder stabil wird.

LG,
Bernhard
---
Programmiere stets so, dass die Maxime Deines Programmierstils Grundlage allgemeiner Gesetzgebung sein könnte
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von EDi »

Finde ich unschön.
Fakt ist dass die beta-Verteilung nicht zu den Daten passt. (und ich bin kein großer Freund von partiellen Transformationen - bigben hat u. A. einen Grund genannt).
(wieso gibt es die Nullen / Einsen?)

Alternativen wären zero/one inflated beta, eine globale transformation https://stats.stackexchange.com/a/134297/1050 oder was ganz anderes.
Vielleicht beschreibst du auch wie der Prozess ist der die Daten produziert (hilf mir oft ein Gefühl zu bekommen wie die Verteilung aussieht, bzw. zu simulieren ).

Ich würde an die Rohdaten anpassen.
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.
Curnen
Beiträge: 27
Registriert: Fr Nov 18, 2016 3:45 pm

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von Curnen »

EDi hat geschrieben: Mi Nov 15, 2017 5:35 pm Alternativen wären zero/one inflated beta, eine globale transformation https://stats.stackexchange.com/a/134297/1050 oder was ganz anderes.
Die im dem StackExchange-Beitrag genannte Formel für die Transformation sieht gut aus.
EDi hat geschrieben: Mi Nov 15, 2017 5:35 pm Finde ich unschön.
Fakt ist dass die beta-Verteilung nicht zu den Daten passt.
[...]
Ich würde an die Rohdaten anpassen.
Leider wäre das Perlen vor die Säue geworfen. Ich bin Molekularbiologe in einem Labor, wo Doktoranden daran gemessen werden, wie viele Proben pro Stunde sie pipettieren können. Dass ich zufällig ein klein wenig "Computerkram" kann ist nett (deswegen mache ich ihn auch nicht nur für mich sondern auch regelmäßig für andere Leute), aber Analysen/Modelle zu verbesseren oder gar neue Programmpakete für Daten zu schreiben gehört nicht in meinen Aufgabenbereich oder zum Gebiet unserer Arbeitsgruppe. Hobbys habe ich schon genug ;-)

Im Übrigen habe ich für meine Auswertung der Daten dann letztlich das GAM-Modell benutzt, das du mir vor drei Jahren mal vorgeschlagen hattest. Hierbei geht es letztlich nur darum für die Dissertation zu zeigen, dass das GAM-Verfahren genauere Ergebnisse liefert das das Standardpaket "Methylseeker" auf Bioconductor. Das kriege ich aber nicht zum Laufen bzw. macht so viele Probleme (u.a. muss man die Daten vorher in ein unübersichtlich verschachteltes Ausgangsformat packen), dass ich mir in der Originalpublikation angeschaut habe, was es tut und das ist letztlich nichts anderes als eine Beta-Funktion zu fitten.
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von EDi »

Im Übrigen habe ich für meine Auswertung der Daten dann letztlich das GAM-Modell benutzt, das du mir vor drei Jahren mal vorgeschlagen hattest. Hierbei geht es letztlich nur darum für die Dissertation zu zeigen, dass das GAM-Verfahren genauere Ergebnisse liefert das das Standardpaket "Methylseeker" auf Bioconductor. Das kriege ich aber nicht zum Laufen bzw. macht so viele Probleme (u.a. muss man die Daten vorher in ein unübersichtlich verschachteltes Ausgangsformat packen), dass ich mir in der Originalpublikation angeschaut habe, was es tut und das ist letztlich nichts anderes als eine Beta-Funktion zu fitten.
Ich erinnere mich wage...
Hab mir mal Methyl-seeker angeschaut. Letzendenlich ist das sowas wie ein Mixture Model zum Klassifizieren (zwei Gruppen mit unterschiedlichen Verteilungen).
Dazu könntest du dir auch mal flexmix::betamix anschauen, um ein mixture model mit 2 beta-Verteilungen zu fitten (https://stats.stackexchange.com/questio ... ll-example). Vermutlich besser als GAM, wenn es dir um das Klassifizieren geht...?
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.
Curnen
Beiträge: 27
Registriert: Fr Nov 18, 2016 3:45 pm

Re: Beta-Funktion auf Methylomdaten fitten

Beitrag von Curnen »

Nachdem ich es ein paar Tage liegen lassen musste, kann ich jetzt vermelden: Es läuft! Vielen vielen Dank nochmal für eure Hilfe!

functions_to_fit.R

Code: Alles auswählen


# Modified Beta function with shape1=shape2 fixed as used with Methylseeker
dbetamod <- function (x, shape1, ncp = 0, log = FALSE) {
  dbeta(x,shape1 = shape1, shape2 = shape1, ncp = ncp, log = log)
}

# Define the Kumaraswamy distribution's PDF
pkumaraswamy <- function(x, a, b) {
  a * b * x ^ (a - 1) * (1 - x ^ a) ^ (b - 1)
}

################################
# Fitting distribution functions
################################


fit_dbeta <- function(x,
                      startvalues=c(0.9,0.6),
                      interval,
                      sample) {
  
  if(any(is.na(x))){
  
      return(NA)
  
    } else {
    
    
    out <- tryCatch({
      message(paste("Fitting Beta-Distribution on interval", interval, "of", sample))
       x <- (x * (length(x)-1) + 0.5) / length(x)
     
      # delete special chars
      sample <- gsub("[^[:alnum:]]","",sample,perl=TRUE)
      
      library("MASS")
      
      fit.beta <-
        fitdistr(x,
                 dbeta,
                 start = list(shape1 = startvalues[1], shape2 = startvalues[2]))
      
      #evaluate fit.beta to return it
      list(fit.beta)
      
    },
    error = function(cond) {
      message(paste("Fitting failed on interval", interval, "of", sample))
      message(cond)
      return(NA)
    })
    return(out)
  }
}

################


fit_dbetamod <- function(x,
                      startvalues=c(0.7),
                      interval,
                      sample) {
  
  if(any(is.na(x))){
    
    return(NA)
    
  } else {
    
    out <- tryCatch({
      message(paste("Fitting restricted Beta-Distribution on interval", interval, "of", sample))
            x <- (x * (length(x)-1) + 0.5) / length(x)
      
      # delete special chars
      sample <- gsub("[^[:alnum:]]","",sample,perl=TRUE)
      
      library("MASS")
      
      fit.betamod <-
        fitdistr(x,
                 dbetamod,
                 start = list(shape1 = startvalues[1]),
                 method = "Brent",
                 lower = .001, 
                 upper = 4)
 
      #evaluate fit.betamod to return it
      list(fit.betamod)
      
    },
    error = function(cond) {
      message(paste("Fitting failed on interval", interval, "of", sample))
      message(cond)
      return(NA)
    })
    return(out)
  }
}

################


fit_kuma <- function(x,
                         startvalues=c(0.9,0.7),
                         interval,
                         sample) {
  if(any(is.na(x))){
    
    return(NA)
    
  } else {
    
    out <- tryCatch({
      message(paste("Fitting Kumaraswamy-Distribution on interval", interval, "of", sample))
      x <- (x * (length(x)-1) + 0.5) / length(x)
      
      # delete special chars
      sample <- gsub("[^[:alnum:]]","",sample,perl=TRUE)
      
      library("MASS")
      
      fit.kuma <-
        fitdistr(x,
                 pkumaraswamy,
                 start = list(a = startvalues[1],b = startvalues[2]))

      #evaluate fit.kuma to return it
      list(fit.kuma)
      
    },
    error = function(cond) {
      message(paste("Fitting failed on interval", interval, "of", sample))
      message(cond)
      # Choose a return value in case of error
      return(NA)
    })
    return(out)
  }
}

Code: Alles auswählen

source("./functions_to_fit.R")

myname <- "SampleA"
methdata <- list(sample(seq(0,1,by=0.0625),1e4, replace = TRUE, prob = c(40,1,2,1,4,1,2,1,8,1,2,1,4,1,2,1,100)),
                 sample(seq(0,1,by=0.0625),1e4, replace = TRUE, prob = c(30,1,2,1,4,1,2,1,30,1,2,1,4,1,2,1,20)),
                 sample(seq(0,1,by=0.0625),1e4, replace = TRUE, prob = c(70,1,2,1,4,1,2,1,3,1,2,1,4,1,2,1,40)))

fits_beta_methdata <- mapply(fit_dbeta,
                             x=methdata,
                             interval=c(1:length(methdata)),
                             sample=rep(myname,length(methdata)))


fits_betamod_methdata <- mapply(fit_dbetamod,
                                x=methdata,
                                interval=c(1:length(methdata)),
                                sample=rep(myname,length(methdata)))


fits_kuma_methdata <- mapply(fit_kuma,
                             x=methdata,
                             interval=c(1:length(methdata)),
                             sample=rep(myname,length(methdata)))

parameters <- list(
  "parameters_beta" = do.call("rbind",lapply(fits_beta_methdata,"[[",1)),
  "parameters_betamod" = do.call("rbind",lapply(fits_betamod_methdata,"[[",1)),
  "parameters_kuma" = do.call("rbind",lapply( fits_kuma_methdata,"[[",1)))

Antworten