Mercurial > repos > jfrancoismartin > mixmodel4repeated_measures

comparison diagmfl.R @ 0:1422de181204 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/mixmodel4repeated_measures commit 6ea32b3182383c19e5333201d2385a61d8da3d50
author jfrancoismartin
date Wed, 10 Oct 2018 05:18:42 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1422de181204
1 ###############################################################################################
2 ## Diagnostics graphics pour les modèles linéaires mixtes de type "mfl" ##
3 ############################################################################################### ##
4 ## Input: ##
5 ## mfl, un modèle linéaire mixte généré par le module lmixedm de JF Martin (fonction ##
6 ## lmRepeated2FF), i.e. : un modèle linéaire mixte de type lmer créé par la formule ##
7 ## mfl <- lmer( vd ~ time + fixfact + time:fixfact + (1| subject), ids) ##
8 ## => noms de colonnes importants ##
9 ## => 1 seul effet aléatoire sur l'intercept ##
10 ###############################################################################################
11 ## Output : ##
12 ## Les graphics, les calculs associés et les notations* utilisées dans le script suivent ##
13 ## l'article de Singer et al (2016) Graphical Tools for detedcting departures from linear ##
14 ## mixed model assumptions and some remedial measures, International Statistical Review ##
15 ## (doi:10.1111/insr.12178) ##
16 ## * ajout d'une ou 2 lettres de typage de la variables ##
17 ##
18 ## Script adapté de http://www.ime.unicamp.br/~cnaber/residdiag_nlme_v22.R pour fonction- ##
19 ## ner avec un modèle lmer (et non lme), des sujets avec des identifiants non numériques, ##
20 ## et des observations non ordonnées sujet par sujet (dernier point à vérifier.) ##
21 ## ############################################################################################
22 ## Remarques sur les calculs numériques ##
23 ## - l'inverse d'une matrice est calculée à partir de la fonction ginv du package MASS ##
24 ## (Moore- Penrose generalized inverse) au lieu de la fonction "solve" ##
25 ## - la racine carrée des matrices sont calculées grâce à une déomposition SVD (car ##
26 ## s'applique normalement ici uniquement à des matrices symétriques positives). A ##
27 ## remplacer par la fonction sqrtm{expm} si des erreurs apparaissent ??? ##
28 ###############################################################################################
29
30
31 library(ggplot2)
32 library(gridExtra)
33 library(grid)
34
35
36 ##-------------------------------------------------------------------------------------------------##
37 ## Helpers
38 ##-------------------------------------------------------------------------------------------------##
39
40 ## square root of a matrix
41 ## http://www.cs.toronto.edu/~jepson/csc420/notes/introSVD.pdf (page 6)
42 ## (for matMN a n x n matrix that symetric and non-negative definite)
43
44 sqrtmF<-function(matMN){
45 matMN <- as.matrix(matMN)
46 # ## check that matMN is symetric
47 # if(!all(t(matMN==matMN)))
48 # stop("matMN must be symetric.")
49 svd_dec <- svd(matMN)
50 invisible(svd_dec$u%*%sqrt(diag(svd_dec$d))%*%t(svd_dec$v))
51 }
52
53
54 ## qqplotF
55 ## adapted from https://gist.github.com/rentrop/d39a8406ad8af2a1066c
56
57
58 qqplotF <- function(x,
59 distribution = "norm", ...,
60 line.estimate = NULL,
61 conf = 0.95,
62 labels = names(x)){
63 q.function <- eval(parse(text = paste0("q", distribution)))
64 d.function <- eval(parse(text = paste0("d", distribution)))
65 x <- na.omit(x)
66 ord <- order(x)
67 n <- length(x)
68 P <- ppoints(length(x))
69 daf <- data.frame(ord.x = x[ord], z = q.function(P, ...))
70
71 if(is.null(line.estimate)){
72 Q.x <- quantile(daf$ord.x, c(0.25, 0.75))
73 Q.z <- q.function(c(0.25, 0.75), ...)
74 b <- diff(Q.x)/diff(Q.z)
75 coef <- c(Q.x[1] - b * Q.z[1], b)
76 } else {
77 coef <- coef(line.estimate(ord.x ~ z))
78 }
79
80 zz <- qnorm(1 - (1 - conf)/2)
81 SE <- (coef[2]/d.function(daf$z,...)) * sqrt(P * (1 - P)/n)
82 fit.value <- coef[1] + coef[2] * daf$z
83 daf$upper <- fit.value + zz * SE
84 daf$lower <- fit.value - zz * SE
85
86 if(!is.null(labels)){
87 daf$label <- ifelse(daf$ord.x > daf$upper | daf$ord.x < daf$lower, labels[ord],"")
88 }
89
90 p <- ggplot(daf, aes(x=z, y=ord.x)) +
91 geom_point() +
92 geom_abline(intercept = coef[1], slope = coef[2], col = "red") +
93 geom_line(aes(x=z, y = lower),daf, col = "red", linetype = "dashed") +
94 geom_line(aes(x=z, y = upper),daf, col = "red", linetype = "dashed") +
95 #geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.2)+
96 xlab("")+ylab("")
97 if(!is.null(labels)) p <- p + geom_text( aes(label = label))
98
99 return(p)
100 #print(p)
101 #coef
102 }
103
104
105 ## histogramm
106 histF <- function(x, sd_x = NULL, breaks = "scott"){
107
108 if(is.null(sd_x)){ sd_x <- sd(x)}
109
110 ## Bandwith estimation (default is Scott)
111 if(!breaks %in% c("sqrt", "sturges", "rice", "scott", "fd"))
112 breaks <- "scott"
113
114 if(breaks %in% c("sqrt", "sturges", "rice")){
115 k <- switch(breaks,
116 sqrt = sqrt(length(x)),
117 sturges = floor(log2(x))+1,
118 rice = floor(2*length(x)^(1/3))
119 )
120 bw <- diff(range(x))/k
121 }else{
122 bw <- switch(breaks,
123 scott = 3.5*sd_x/length(x)^(1/3),
124 fd = diff(range(x))/(2*IQR(x)/length(x)^(1/3))
125 )
126 }
127
128
129 daf <- data.frame(x=x)
130 ## graph
131 return(ggplot(data=daf, aes(x)) +
132 geom_histogram(aes(y = ..density..),
133 col="black", fill = "grey", binwidth = bw)+
134 geom_density(size = 1.2,
135 col = "blue",
136 linetype = "blank",
137 fill = rgb(0,0,1, 0.1))+
138 stat_function(fun = dnorm,
139 args = list(mean = 0, sd = sd_x),
140 col = "blue", size = 1.2)+
141 theme(legend.position="none")+
142 xlab(""))
143 }
144
145
146
147 ##-------------------------------------------------------------------------------------------------##
148 ## Main function
149 ##-------------------------------------------------------------------------------------------------##
150
151
152 diagmflF <- function(mfl, title = "",
153 outlier.limit = 3,
154 pvalCutof = 0.05,
155 least.confounded = FALSE){
156
157 ## initialisation -------------------------------------------------------------------------------------
158 responseC<- "vd"
159 unitC <- "subject"
160 timeC<- "time"
161 fixfactC<-"fixfact"
162 hlimitN <- outlier.limit
163
164
165
166 ## extracting information from mfl models -------------------------------------------------------------
167 df <- mfl@frame
168 yVn <- df[, responseC]
169 nobsN <- length(yVn)
170 idunitVc <- unique(df[, unitC])
171 nunitN <- length(unique(idunitVc))
172 xMN <- mfl@pp$X
173 pN <- ncol(xMN)
174 zMN <- as.matrix(t(mfl@pp$Zt))
175 gMN <- as.matrix(as.data.frame(VarCorr(mfl))[1, "vcov"]) ## Valable seulement pour 1 seul effet aléatoire
176 gammaMN <- as.matrix(kronecker(diag(nunitN), gMN))
177 sigsqN <- as.data.frame(VarCorr(mfl))[2, "vcov"]
178 rMN <- sigsqN*diag(nobsN)
179 vMN <- (zMN%*%gammaMN%*%t(zMN)) + rMN
180 invvMN <- MASS::ginv(vMN)
181 hMN <- MASS::ginv(t(xMN)%*%invvMN%*%xMN)
182 qMN <- invvMN - invvMN%*%xMN%*%(hMN)%*%t(xMN)%*%invvMN
183 eblueVn<-mfl@beta
184 eblupVn<-gammaMN%*%t(zMN)%*%invvMN%*%(yVn-xMN%*%eblueVn) ## equivalent de ranef(mfl)
185 rownames(eblupVn) <- colnames(zMN)
186
187
188 ## Calculs of matrices and vectors used in graph diagnosics ---------------------------------------------
189
190 ## Marginal and individual predictions, residuals and variances
191 marpredVn <- xMN%*%eblueVn
192 marresVn <- yVn - marpredVn
193 marvarMN <- vMN - xMN%*%hMN%*%t(xMN)
194
195 condpredVn <- marpredVn + zMN%*%eblupVn
196 condresVn <- yVn - condpredVn
197 condvarMN <- rMN%*%qMN%*%rMN
198
199
200 ## Analysis of marginal and conditional residuals
201 stmarresVn <-stcondresVn <- rep(0,nobsN)
202 lesverVn <- rep(0, nunitN)
203 names(lesverVn )<- idunitVc
204
205 for(i in 1:nunitN){
206
207 idxiVn <- which(df[, unitC] == idunitVc[i]) ## position des observations du sujet i
208 miN <- length(idxiVn)
209
210 ## standardization of marginal residual
211 stmarresVn[idxiVn] <- as.vector(solve(sqrtmF(marvarMN[idxiVn,idxiVn]))%*%marresVn[idxiVn])
212
213 ##Standardized Lessafre and Verbeke's measure
214 auxMN <- diag(1, ncol = miN, nrow =miN)- stmarresVn[idxiVn]%*%t(stmarresVn[idxiVn])
215 lesverVn[i] <- sum(diag(auxMN%*%t(auxMN)))
216
217 ## standardization of conditional residual
218 stcondresVn[idxiVn] <- as.vector(solve(sqrtmF(condvarMN[idxiVn,idxiVn]))%*%condresVn[idxiVn])
219 }
220 lesverVn <- lesverVn/sum(lesverVn)
221
222
223 ## Least confounded residuals (à valider !)
224 ## résultats différents des valeurs estimées via le script de Singer, car
225 ## non unicité des vecteurs propres de la décomposition spectrale ?
226 if(least.confounded){
227 sqrMN <- sqrtmF(rMN)
228 specDec <- eigen((sqrMN%*%qMN%*%sqrMN), symmetric = TRUE, only.values = FALSE)
229
230 cMN <- t(sqrt(solve(diag((specDec$values[1:(nobsN -pN)])))) %*% t(specDec$vectors[1:(nobsN -pN),1:(nobsN -pN)]) %*%
231 solve(sqrtmF(rMN[1:(nobsN -pN),1:(nobsN -pN)])) )
232
233 lccondresVn <- (cMN%*%condresVn[1:(nobsN -pN)])/sqrt(diag(cMN%*%condvarMN[1:(nobsN -pN),1:(nobsN -pN)]%*%t(cMN)))
234 }
235
236
237 ## EBLUP analysis (Mahalanobis' distance)
238 varbMN <- gammaMN%*%t(zMN)%*%qMN%*%zMN%*%gammaMN
239 mdistVn <- rep(0, nunitN)
240 for(i in 1:nunitN){
241 mdistVn[i] <- eblupVn[i]^2/varbMN[i, i]
242 }
243 mdistVn <- mdistVn/sum(mdistVn)
244
245
246 ## Combine data and results in 2 data frames for easy plotting with ggplot2 -----------------------------
247
248 ## long data frame (all observations)
249
250 df <- data.frame(df,
251 mar.pred = marpredVn,
252 mar.res = marresVn,
253 st.mar.res = stmarresVn,
254 cond.pred = condpredVn,
255 cond.res = condresVn,
256 st.cond.res = stcondresVn)
257
258 if(!("fixfact" %in% colnames(df)))
259 df$fixfact <- rep(1, nrow(df))
260 df$numTime <- as.numeric(levels(as.factor(df$time)))
261
262
263 ## short data frame (1 row per unit)
264
265 unitDf <- data.frame(unit = idunitVc,
266 eblup = eblupVn,
267 lvm = lesverVn,
268 mal = mdistVn)
269
270 unitDf$fixfact <- sapply(1:nrow(unitDf),
271 function(i){
272 unique(df[which(df[, unitC] == unitDf$unit[i]),
273 fixfactC])
274 })
275
276 unitDf$se <- rep(NA, nrow(unitDf))
277 re <- ranef(mfl, condVar =TRUE)
278 for(i in 1:nrow(unitDf))
279 unitDf$se[i] <- sqrt(attr(re[[unitC]], "postVar")[1,1,i])
280 unitDf$upper <- unitDf$eblup+unitDf$se*1.96
281 unitDf$lower <- unitDf$eblup-unitDf$se*1.96
282
283
284 ## Outliers "annotations"
285 df$marres.out <- rep("", nrow(df))
286 df$marres.out[abs(df$st.mar.res)>hlimitN] <-
287 paste(df[abs(df$st.mar.res)>hlimitN, unitC],
288 df[abs(df$st.mar.res)>hlimitN, timeC],
289 sep = ".")
290 df$marres.out <- paste(" ", df$marres.out, sep ="")
291
292 df$condres.out <- rep("", nrow(df))
293 df$condres.out[abs(df$st.cond.res)>hlimitN] <-
294 paste(df[abs(df$st.cond.res)>hlimitN, unitC],
295 df[abs(df$st.cond.res)>hlimitN, timeC],
296 sep = ".")
297 df$condres.out <- paste(" ", df$condres.out, sep ="")
298
299 ## Diagnostic Plots -------------------------------------------------------------------------------------
300
301 ## Linearity of effect and outlying observations
302 p1 <- ggplot(data = df,
303 aes(x=mar.pred,
304 y=st.mar.res,
305 colour=fixfact)) +
306 geom_point(size =2) +
307 geom_hline(yintercept = 0, col = "grey")+
308 geom_smooth(aes(x=mar.pred, y=st.mar.res), data = df, se = FALSE, col = "blue", method = "loess")+
309 ggtitle("Linearity of effects/outlying obervations")+
310 xlab("Marginal predictions")+
311 ylab("Standardized marginal residuals")+
312 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+
313 geom_hline(yintercept = c(-1,1)*hlimitN, linetype = "dashed")+
314 geom_text(aes(label = marres.out, col = fixfact), hjust=0, vjust=0)
315
316 # p1hist <- histF(df$st.mar.res, sd_x =1)+
317 # xlab("Standardized marginal residuals")
318
319
320 ## Presence of outlying observations and homoscedacity of residuals
321 p2 <- ggplot(data = df,
322 aes(x=cond.pred,
323 y=st.cond.res,
324 colour=fixfact)) +
325 geom_point(size =2) +
326 geom_hline(yintercept = 0, col = "grey")+
327 geom_smooth(aes(x=cond.pred, y=st.cond.res), data = df, se = FALSE, col = "blue", method = "loess")+
328 ggtitle("Homosedasticity of conditional residuals/outlying observations")+
329 xlab("Individual predictions")+
330 ylab("Standardized conditional residuals")+
331 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+
332 geom_hline(yintercept = c(-1,1)*hlimitN, linetype = "dashed")+
333 geom_text(aes(label = condres.out, col = fixfact), hjust=0, vjust=0)
334
335 # p2hist <- histF(df$st.cond.res, sd_x =1)+
336 # xlab("Standardized conditional residuals")
337
338
339 ## Normality of residuals
340 if(least.confounded){
341 p3 <- qqplotF(x = lccondresVn,
342 distribution = "norm",
343 line.estimate = NULL,
344 conf = 0.95)+
345 xlab("Standard normal quantiles")+
346 ylab("Least confounded conditional residual quantiles")+
347 ggtitle("Normality of conditional error (least confounded)")+
348 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))
349
350 # p3hist <- histF(lccondresVn, sd_x =1)+
351 # xlab("Least confounded conditional residuals")
352 }else{
353 p3 <-qqplotF(x = df$st.cond.res,
354 distribution = "norm",
355 line.estimate = NULL,
356 conf = 0.95)+
357 xlab("Standard normal quantiles")+
358 ylab("Standardized conditional residual quantiles")+
359 ggtitle("Normality of conditional error")+
360 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))
361 # p3hist <-p2hist
362 }
363
364 ## Within-units covariance structure
365 p4 <- ggplot(data = unitDf,
366 aes(x=unit,
367 y=lvm,
368 colour=fixfact)) +
369 geom_point(size =2) +
370 theme(legend.position="none")+
371 xlab(unitC)+
372 ylab("Standardized Lesaffre-Verbeke measure")+
373 geom_hline(yintercept = 2*mean(unitDf$lvm), linetype = "dashed")+
374 geom_text(aes(label = unit, col = fixfact),
375 data = unitDf[unitDf$lvm>2*mean(unitDf$lvm), ],
376 hjust=0, vjust=0)+
377 ggtitle("Within-units covariance matrice")+
378 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))
379
380 ## EBLUP modif lmerTest v3 mais plante si pas d'effet aléatoire
381 #pvl <- ranova(model=mfl, reduce.terms = TRUE)
382 #pvrnd <- pvl[[6]][2]
383 #ggtitle(paste("Random effect on intercept (LRT p-value = ",round(pvrnd,digits = 5), ")", sep = ""))+
384
385 p5 <-
386 ggplot(aes(x = eblup, y = unit, col = fixfact), data = unitDf)+
387 geom_point(size =3)+
388 geom_segment(aes(xend = lower, yend = unit)) +
389 geom_segment(aes(xend = upper, yend = unit))+
390 ggtitle("Random effect on intercept")+
391 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+
392 geom_vline(xintercept = 0, linetype = "dashed")+
393 ylab(unitC)
394
395 ## p6
396 p6 <-qqplotF(x = unitDf$mal,
397 distribution = "chisq",
398 df= 1,
399 line.estimate = NULL,
400 conf = 0.95)+
401 xlab("Chi-squared quantiles")+
402 ylab("Standadized Mahalanobis distance")+
403 ggtitle("Normality of random effect")+
404 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))
405
406 ## Outlying subjects
407 p7 <-
408 ggplot(aes(y = mal, x= unit, col = fixfact), data = unitDf)+
409 geom_point(size =3)+
410 ylab("Standardized Mahalanobis distance")+
411 geom_vline(xintercept = 0, linetype = "dashed")+
412 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+
413 geom_hline(yintercept = 2*mean(unitDf$mal), linetype = "dashed")+
414 geom_text(aes(label = unit, col = fixfact),
415 data = unitDf[unitDf$mal>2*mean(unitDf$mal), ],
416 hjust=1, vjust=0)+
417 ggtitle("Outlying subjects")+
418 xlab(unitC)
419
420 ## "Data" and "modeling" Plots --------------------------------------------------------------------------
421
422 ## Individual time-course
423 rawPlot <- ggplot(data = df,
424 aes(x=time, y=vd, colour=fixfact, group=subject)) +
425 geom_line() + ggtitle("Individual time-courses ")+
426 theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))
427
428 ## Post-hoc estimates (modification due to lmerTest v3)
429 ddlsm1 <- data.frame(difflsmeans(mfl,test.effs=NULL))
430 colnames(ddlsm1)[9] <- "pvalue"
431 # ddlsm1$name <- sapply(rownames(ddlsm1),
432 # function(nam){
433 # strsplit(nam, split = " ", fixed =TRUE)[[1]][1]
434 # })
435 # ddlsm1$detail <- sapply(rownames(ddlsm1),
436 # function(nam){
437 # paste(strsplit(nam, split = " ", fixed =TRUE)[[1]][-1],
438 # collapse= "")
439 # })
440 #
441 # colnames(ddlsm1)<- make.names(colnames(ddlsm1))
442 ddlsm1$Significance <- rep("NS", nrow(ddlsm1))
443
444 ## modif JF pour tenir compte du seuil de pvalues defini par le user
445 ddlsm1$Significance[which(ddlsm1$pvalue <pvalCutof & ddlsm1$pvalue >=0.01)] <- "p-value < threshold"
446 ddlsm1$Significance[which(ddlsm1$pvalue <0.01 & ddlsm1$pvalue >=0.005)] <- "p-value < 0.01"
447 ddlsm1$Significance[which(ddlsm1$pvalue <0.005)] <- "p-value < 0.005"
448
449 phPlot <- ggplot(ddlsm1, aes(x = levels, y = Estimate))+
450 facet_grid(facets = ~term, ddlsm1,scales = "free", space = "free")+
451 geom_bar( aes(fill = Significance), stat="identity")+
452 theme(axis.text.x = element_text(angle = 90, hjust = 1))+
453 scale_fill_manual(
454 values = c("NS" = "grey",
455 "p-value < threshold" = "yellow",
456 "p-value < 0.01" = "orange",
457 "p-value < 0.005" = "red"))+
458 geom_errorbar(aes(ymin = lower, ymax =upper ), width=0.25)+
459 ggtitle("Post-hoc estimates")+xlab("")+
460 theme(plot.title = element_text(size = rel(1.2), face = "bold"))
461
462 ## Final plotting
463
464 blank<-rectGrob(gp=gpar(col="white"))
465 rectspacer<-rectGrob(height = unit(0.1, "npc"), gp=gpar(col="grey"))
466
467 grid.arrange(blank,
468 rawPlot, blank, phPlot,
469 rectspacer,
470 p1,blank, p2, blank, p3, blank, p4,
471 blank,
472 p5,blank, p6, blank, p7, blank,
473 blank,blank,
474 top = textGrob(title,gp=gpar(fontsize=40,font=4)),
475 layout_matrix = matrix(c(rep(1,7),
476 2, 3, rep(4,3), 20,21,
477 rep(5,7),
478 6:12,
479 rep(13,7),
480 14:18, rep(19,2)),
481 ncol=7, nrow=6, byrow=TRUE),
482 heights= c(0.1/3, 0.3, 0.1/3, 0.3, 0.1/3, 0.3 ),
483 widths = c(0.22, 0.04, 0.22,0.04 , 0.22, 0.04, 0.22))
484
485 invisible(NULL)
486
487 }
488
489 #diagmflF(mfl, title = "", outlier.limit = 3, least.confounded =TRUE)
490
491 plot.res.Lmixed <- function(mfl, df, title = "", pvalCutof = 0.05) {
492
493 nameVar <- colnames(df)[4]
494 fflab <- colnames(df)[1]
495 ## Individual time-course
496 rawPlot <-
497 ggplot(data = df, aes(x=df[[2]], y=df[[4]], colour=df[[1]], group=df[[3]])) +
498 geom_line() + ggtitle("Individual time-courses (raw data)")+
499 ylab(nameVar) +
500 xlab(label = colnames(df)[2])+
501 theme(legend.title = element_blank() , legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))
502
503 ## Boxplot of fixed factor
504 bPlot <-
505 ggplot(data = df, aes(y=df[[4]], x=df[[1]], color=df[[1]]))+
506 geom_boxplot(outlier.colour="red", outlier.shape=8, outlier.size=4)+
507 ggtitle(paste("Boxplot by ",fflab,sep=""))+xlab("")+ylab("")+
508 theme(legend.title = element_blank(), plot.title = element_text(size = rel(1.2), face = "bold"))
509
510 ## Post-hoc estimates
511
512 ddlsm1 <- mfl
513 ddlsm1$name <- rownames(ddlsm1)
514 # ddlsm1$name <- sapply(rownames(ddlsm1),
515 # function(nam){
516 # strsplit(nam, split = " ", fixed =TRUE)[[1]][1]
517 # })
518 # ddlsm1$detail <- sapply(rownames(ddlsm1),
519 # function(nam){
520 # paste(strsplit(nam, split = " ", fixed =TRUE)[[1]][-1],
521 # collapse= "")
522 # })
523 #
524 #colnames(ddlsm1)<- make.names(colnames(ddlsm1))
525 ddlsm1$Significance <- rep("NS", nrow(ddlsm1))
526
527 ## modif JF pour tenir compte du seuil de pvalues defini par le user
528 options("scipen"=100, "digits"=5)
529 pvalCutof <- as.numeric(pvalCutof)
530 bs=0.05; bm=0.01; bi=0.005
531 if (pvalCutof >bm) {bs <- pvalCutof} else
532 if (pvalCutof <bm & pvalCutof >bi) {bm <- pvalCutof; bs <- pvalCutof} else
533 if (pvalCutof < bi) {bi <- pvalCutof; bm <- pvalCutof; bs <- pvalCutof}
534 lbs <- paste("p-value < ",bs, sep="")
535 lbm <- paste("p-value < ",bm, sep="")
536 lbi <- paste("p-value < ",bi, sep="")
537
538 cols <- paste("p-value < ",bs, sep="")
539 colm <- paste("p-value < ",bm, sep="")
540 coli <- paste("p-value < ",bi, sep="")
541 valcol <- c("grey","yellow","orange","red")
542 names (valcol) <- c("NS",lbs,lbm,lbi)
543 ddlsm1$Significance[which(ddlsm1$p.value<= bs)] <- lbs
544 ddlsm1$Significance[which(ddlsm1$p.value<bs & ddlsm1$p.value>= bm)] <- lbm
545 ddlsm1$Significance[which(ddlsm1$p.value<bi)] <- lbi
546
547 phPlot <-
548 ggplot(ddlsm1, aes(x = levels, y = Estimate))+
549 facet_grid(facets = ~term, ddlsm1,scales = "free", space = "free")+
550 geom_bar( aes(fill = Significance), stat="identity")+
551 theme(axis.text.x = element_text(angle = 90, hjust = 1))+
552 scale_fill_manual(
553 # values = c("NS" = "grey", "p-value < threshold" = "yellow","p-value < 0.01" = "orange","p-value < 0.005" = "red"))+
554 # values = c("NS" = 'grey', "pvalue < 0.05 "= 'yellow',"p-value < 0.01" = 'orange',"p-value < 0.005" = 'red'))+
555 # values = c("NS = grey", "p-value < threshold = yellow","p-value < 0.01 = orange","p-value < 0.005 = red"))+
556 values = valcol )+
557 geom_errorbar(aes(ymin = Lower.CI, ymax =Upper.CI ), width=0.25)+
558 # ggtitle(paste("Post-hoc estimates with p-value threshold = ",pvalCutof,sep=""))+xlab("")+
559 ggtitle("Post-hoc estimates ")+xlab("")+
560 theme(plot.title = element_text(size = rel(1.2), face = "bold"))
561
562 ## Final plotting
563 grid.arrange(arrangeGrob(rawPlot,bPlot,ncol=2),
564 phPlot,nrow=2,
565 top = textGrob(title,gp=gpar(fontsize=32,font=4))
566 )
567
568 }