Mercurial > repos > marie-tremblay-metatoul > asca

comparison ssq.R @ 0:93312041f1d5 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/ascaw4m commit 7ea9b0f8abc5a60c2c04fd2098788497f14766b6
author marie-tremblay-metatoul
date Fri, 21 Sep 2018 05:51:14 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:93312041f1d5
1 library(lmdme)
2 library(MetStaT)
3
4 ## Data : attention standardiser au prealable
5 data <- read.table("E:/PROJETS/Asca_W4M/Test_Matlab_data.txt", sep="\t", dec=",", header=TRUE, row.names=1)
6 design <- read.table("E:/PROJETS/Asca_W4M/Test_Matlab_design.txt", sep="\t", dec=",", header=TRUE)
7 design[,1] <- as.factor(design[,1])
8 design[,2] <- as.factor(design[,2])
9
10 ## Verifier noms
11
12
13 fit <- lmdme(model=~F1 + F2 + F1:F2, data=data, design=design)
14
15 permuted <- permutation(model=~F1*F2, data=data, design=design, NPermutations=100, nCpus=3)
16
17 decomposition(fit, decomposition = "pca", scale="none", type="coefficient")
18
19
20
21
22 ssq <- function(fit)
23 {
24 Overall_means <- sum(sum(fitted.values(fit)$'(Intercept)'^2))/sum(sum(data^2))
25 Factors <- c(sum(sum(fitted.values(fit)$'F1'^2))/sum(sum(data^2)), sum(sum(fitted.values(fit)$'F2'^2))/sum(sum(data^2)))
26 Interactions <- sum(sum(fitted.values(fit)$'F1:F2'^2))/sum(sum(data^2))
27 Residuals <- 1 - Overall_means - Factors[1] - Factors[2] - Interactions
28
29 return(list(Overall_means, Factors, Interactions, Residuals))
30 }
31
32
33 par(mfrow=c(2,2))
34 biplot(fit, xlabs="o", mfcol=NULL)
35 ##Just the term of interest
36 biplot(fit, xlabs="o", term="F1")
37 ##In separate graphics
38 biplot(fit, xlabs="o", term=c("F1", "F2"), mfcol=c(1,1))
39 ##All terms in the same graphic
40 biplot(fit, xlabs="o", mfcol=c(1,3))
41
42 test <-lapply(permuted, FUN = ssq)
43 test1 <- matrix(unlist(test), ncol=5, byrow=TRUE)
44 test[2:101]
45
46 apply(apply(test1, 2, FUN = function(x){x > x[1]})[-1,], 2 , sum) / (length(test)-1)
47
48 score_moyen <- data.frame(fit@components$F1$rotation)
49
50 score <- data.frame(cbind(design, t(fit@residuals$'F1:F2')%*%fit@components$F1$rotation))
51
52 pc <- fit@components$F1$sdev / sum(fit@components$F1$sdev)
53
54 sp <- ggplot(score_moyen, aes(x=PC1, y=PC2))
55 sp + geom_point(size=2) + xlab(paste("PC1", round(pc[1]*100,1), "%")) + ylab(paste("PC2", round(pc[2]*100,1), "%")) +
56 geom_text(data=score, aes(PC1, PC2, label=Ind, col=F1), size=4, hjust=0, nudge_x=0.05, vjust=0, nudge_y=0.5)
57
58