lme() als aov() Ersatz
Verfasst: Mo Nov 26, 2018 12:20 pm
Ich möchte eine quantitative Variable auf Assoziation mit mehreren kategorischen Variablen testen. Meine erste Idee (Anova) ging leider nicht, weil die Gruppen nicht gleich groß waren. Im Hilfetext zu ?aov() habe ich dann den Hinweis gefunden, man solle für „unbalanced groups“ die Funktion lme() benutzen.
Leider tritt dabei (wieder einmal) das Problem auf, dass ich damit Werkzeuge verwenden muss, deren Komplexität mein Statistikwissen um ein Vielfaches übersteigt. Ich fürchte, ich brauche Hilfe bei der Spezifikation der Modellformel und auch der Interpretation der Ergebnisse.
Mein Datensatz ist im Prinzip so aufgebaut:
Ich habe in 10 Proben (sample_ID) die Expression (sample_expression) von Genen (transcript_ID) bestimmt. Jedes Sample gehört einer Genotype-Gruppe an und jedes Gen fällt in eine Peak Kategorie (peak_status). Jede Sample ID hat also einen festen Genotype und jede Transcript ID einen festen Peak Status.
Wie muss denn meine Modellformel und der PostHoc-Test ausehen, wenn ich auf Assozation mit genotype und peak_status testen will?
Besonders bei dem random-Attribut habe ich bei CrossValidated zum Teil faszinierende Dinge gesehen, wobei leider nicht erklärt wurde, warum es so gemacht wurde und was das bedeuten soll.
Leider beinflusst der Random-Teil die entsprechenden Modelle und damit auch die mögliche statistische Signifikanz deutlich. Basierend auf den Medianen gehe ich davon aus, dass es eine Assoziation mit dem peak_status = "H3K4me3 buffer domain" geben sollte, alles andere sollte vermutlich vernachlässigbar sein. Die echten Daten sehen nämlich so aus:
Ich würde mich sehr freuen, wenn ihr mir da wieder einmal unter die Arme greifen könntet. Vielen Dank!
Matthias
Leider tritt dabei (wieder einmal) das Problem auf, dass ich damit Werkzeuge verwenden muss, deren Komplexität mein Statistikwissen um ein Vielfaches übersteigt. Ich fürchte, ich brauche Hilfe bei der Spezifikation der Modellformel und auch der Interpretation der Ergebnisse.
Mein Datensatz ist im Prinzip so aufgebaut:
Code: Alles auswählen
mydata <- data.frame(
transcript_ID=as.factor(paste("transcript",rep(c(1:50),each=10))),
peak_status=as.factor(rep(c("none","Regular H3K4me3 peak","H3K4me3 buffer domain"),times=c(160,270,70))), # echtes Verhältnis 106:315:1
sample_ID=as.factor(rep(paste("sample",LETTERS[1:10]),times=50)),
genotype=as.factor(rep(rep(c("wt","cchip"),each=5),times=50)),
sample_expression=unlist(lapply(mapply(c,as.list(c(rnorm(43,15,3),rnorm(7,100,20))),as.list(rnorm(50,6,2)),SIMPLIFY=FALSE),function(x){rnorm(10,x[1],x[2])}))
)
Wie muss denn meine Modellformel und der PostHoc-Test ausehen, wenn ich auf Assozation mit genotype und peak_status testen will?
Code: Alles auswählen
library("nlme")
library("multcomp")
fit <- lme(sample_expression ~ genotype + peak_status, data=mydata,random=list(transcript_ID= ~1,sample_ID= ~ 1))
anova(fit)
summary(glht(fit,linfct=mcp(peak_status="Tukey")))
Code: Alles auswählen
# einer von vielen Google-Funden, wo im Random-Teil kuriose Dinge gemacht wurden...
fit <- lme(sample_expression ~ genotype+peak_status, data=mydata,random = ~1|sample_ID/transcript_ID)
Code: Alles auswählen
aggregate(sample_expression ~ genotype + peak_status, data=plotdata_buffered_genes, FUN=median)
genotype peak_status sample_expression
1 wt none 6.388254
2 chip none 5.344367
3 wt Regular H3K4me3 peak 12.299607
4 chip Regular H3K4me3 peak 12.488074
5 wt H3K4me3 buffer domain 147.485491
6 chip H3K4me3 buffer domain 162.046081
anova(fit)
numDF denDF F-value p-value
(Intercept) 1 85542 408.9787 <.0001
genotype 1 85542 12.5218 4e-04
peak_status 2 85542 15.9353 <.0001
summary(glht(fit,linfct=mcp(peak_status="Tukey",genotype = "Tukey")))
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
peak_status: Regular H3K4me3 peak - none == 0 -1.9040 1.9911 -0.956 0.716415
peak_status: H3K4me3 buffer domain - none == 0 75.0342 13.6525 5.496 < 1e-04 ***
peak_status: H3K4me3 buffer domain - Regular H3K4me3 peak == 0 76.9382 13.6775 5.625 < 1e-04 ***
genotype: Dnmt1 -/chip - Dnmt1 +/+ == 0 2.1140 0.5666 3.731 0.000582 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Matthias