lme() als aov() Ersatz

Allgemeine Statistik mit R, die Test-Methode ist noch nicht bekannt, ich habe noch keinen Plan!

Moderatoren: EDi, jogo

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

lme() als aov() Ersatz

Beitrag von Curnen »

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:

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])}))
)
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?

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")))
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.

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)
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:

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)
Ich würde mich sehr freuen, wenn ihr mir da wieder einmal unter die Arme greifen könntet. Vielen Dank!
Matthias
Benutzeravatar
EDi
Beiträge: 1599
Registriert: Sa Okt 08, 2016 3:39 pm

Re: lme() als aov() Ersatz

Beitrag von EDi »

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.

# 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)
Das ist ein sogenannter nested random Effekt und könnte hier sinn machen (mehrere transkripts gehört zu einem sample). Kann man sich so vorstellen, dass zuerst der effekt von sample gemacht wird und wenn man den hat, kommt auf sen wert noch das transkript hinzu. Das kann mam auch als 1|sample + sample:transkript schreiben.

Ich empfehle das glmm buch von Zuur.
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.
Antworten