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?