Mercurial > repos > jfrancoismartin > mixmodel4repeated_measures
diff diagmfl.R @ 1:a3147e3d66e2 draft default tip
Uploaded
| author | melpetera |
|---|---|
| date | Mon, 16 May 2022 12:31:58 +0000 |
| parents | 1422de181204 |
| children |
line wrap: on
line diff
--- a/diagmfl.R Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,568 +0,0 @@ -############################################################################################### -## Diagnostics graphics pour les modèles linéaires mixtes de type "mfl" ## -############################################################################################### ## -## Input: ## -## mfl, un modèle linéaire mixte généré par le module lmixedm de JF Martin (fonction ## -## lmRepeated2FF), i.e. : un modèle linéaire mixte de type lmer créé par la formule ## -## mfl <- lmer( vd ~ time + fixfact + time:fixfact + (1| subject), ids) ## -## => noms de colonnes importants ## -## => 1 seul effet aléatoire sur l'intercept ## -############################################################################################### -## Output : ## -## Les graphics, les calculs associés et les notations* utilisées dans le script suivent ## -## l'article de Singer et al (2016) Graphical Tools for detedcting departures from linear ## -## mixed model assumptions and some remedial measures, International Statistical Review ## -## (doi:10.1111/insr.12178) ## -## * ajout d'une ou 2 lettres de typage de la variables ## -## -## Script adapté de http://www.ime.unicamp.br/~cnaber/residdiag_nlme_v22.R pour fonction- ## -## ner avec un modèle lmer (et non lme), des sujets avec des identifiants non numériques, ## -## et des observations non ordonnées sujet par sujet (dernier point à vérifier.) ## -## ############################################################################################ -## Remarques sur les calculs numériques ## -## - l'inverse d'une matrice est calculée à partir de la fonction ginv du package MASS ## -## (Moore- Penrose generalized inverse) au lieu de la fonction "solve" ## -## - la racine carrée des matrices sont calculées grâce à une déomposition SVD (car ## -## s'applique normalement ici uniquement à des matrices symétriques positives). A ## -## remplacer par la fonction sqrtm{expm} si des erreurs apparaissent ??? ## -############################################################################################### - - -library(ggplot2) -library(gridExtra) -library(grid) - - -##-------------------------------------------------------------------------------------------------## -## Helpers -##-------------------------------------------------------------------------------------------------## - -## square root of a matrix -## http://www.cs.toronto.edu/~jepson/csc420/notes/introSVD.pdf (page 6) -## (for matMN a n x n matrix that symetric and non-negative definite) - -sqrtmF<-function(matMN){ - matMN <- as.matrix(matMN) - # ## check that matMN is symetric - # if(!all(t(matMN==matMN))) - # stop("matMN must be symetric.") - svd_dec <- svd(matMN) - invisible(svd_dec$u%*%sqrt(diag(svd_dec$d))%*%t(svd_dec$v)) -} - - -## qqplotF -## adapted from https://gist.github.com/rentrop/d39a8406ad8af2a1066c - - -qqplotF <- function(x, - distribution = "norm", ..., - line.estimate = NULL, - conf = 0.95, - labels = names(x)){ - q.function <- eval(parse(text = paste0("q", distribution))) - d.function <- eval(parse(text = paste0("d", distribution))) - x <- na.omit(x) - ord <- order(x) - n <- length(x) - P <- ppoints(length(x)) - daf <- data.frame(ord.x = x[ord], z = q.function(P, ...)) - - if(is.null(line.estimate)){ - Q.x <- quantile(daf$ord.x, c(0.25, 0.75)) - Q.z <- q.function(c(0.25, 0.75), ...) - b <- diff(Q.x)/diff(Q.z) - coef <- c(Q.x[1] - b * Q.z[1], b) - } else { - coef <- coef(line.estimate(ord.x ~ z)) - } - - zz <- qnorm(1 - (1 - conf)/2) - SE <- (coef[2]/d.function(daf$z,...)) * sqrt(P * (1 - P)/n) - fit.value <- coef[1] + coef[2] * daf$z - daf$upper <- fit.value + zz * SE - daf$lower <- fit.value - zz * SE - - if(!is.null(labels)){ - daf$label <- ifelse(daf$ord.x > daf$upper | daf$ord.x < daf$lower, labels[ord],"") - } - - p <- ggplot(daf, aes(x=z, y=ord.x)) + - geom_point() + - geom_abline(intercept = coef[1], slope = coef[2], col = "red") + - geom_line(aes(x=z, y = lower),daf, col = "red", linetype = "dashed") + - geom_line(aes(x=z, y = upper),daf, col = "red", linetype = "dashed") + - #geom_ribbon(aes(ymin = lower, ymax = upper), alpha=0.2)+ - xlab("")+ylab("") - if(!is.null(labels)) p <- p + geom_text( aes(label = label)) - - return(p) - #print(p) - #coef -} - - -## histogramm -histF <- function(x, sd_x = NULL, breaks = "scott"){ - - if(is.null(sd_x)){ sd_x <- sd(x)} - - ## Bandwith estimation (default is Scott) - if(!breaks %in% c("sqrt", "sturges", "rice", "scott", "fd")) - breaks <- "scott" - - if(breaks %in% c("sqrt", "sturges", "rice")){ - k <- switch(breaks, - sqrt = sqrt(length(x)), - sturges = floor(log2(x))+1, - rice = floor(2*length(x)^(1/3)) - ) - bw <- diff(range(x))/k - }else{ - bw <- switch(breaks, - scott = 3.5*sd_x/length(x)^(1/3), - fd = diff(range(x))/(2*IQR(x)/length(x)^(1/3)) - ) - } - - - daf <- data.frame(x=x) - ## graph - return(ggplot(data=daf, aes(x)) + - geom_histogram(aes(y = ..density..), - col="black", fill = "grey", binwidth = bw)+ - geom_density(size = 1.2, - col = "blue", - linetype = "blank", - fill = rgb(0,0,1, 0.1))+ - stat_function(fun = dnorm, - args = list(mean = 0, sd = sd_x), - col = "blue", size = 1.2)+ - theme(legend.position="none")+ - xlab("")) -} - - - -##-------------------------------------------------------------------------------------------------## -## Main function -##-------------------------------------------------------------------------------------------------## - - -diagmflF <- function(mfl, title = "", - outlier.limit = 3, - pvalCutof = 0.05, - least.confounded = FALSE){ - - ## initialisation ------------------------------------------------------------------------------------- - responseC<- "vd" - unitC <- "subject" - timeC<- "time" - fixfactC<-"fixfact" - hlimitN <- outlier.limit - - - - ## extracting information from mfl models ------------------------------------------------------------- - df <- mfl@frame - yVn <- df[, responseC] - nobsN <- length(yVn) - idunitVc <- unique(df[, unitC]) - nunitN <- length(unique(idunitVc)) - xMN <- mfl@pp$X - pN <- ncol(xMN) - zMN <- as.matrix(t(mfl@pp$Zt)) - gMN <- as.matrix(as.data.frame(VarCorr(mfl))[1, "vcov"]) ## Valable seulement pour 1 seul effet aléatoire - gammaMN <- as.matrix(kronecker(diag(nunitN), gMN)) - sigsqN <- as.data.frame(VarCorr(mfl))[2, "vcov"] - rMN <- sigsqN*diag(nobsN) - vMN <- (zMN%*%gammaMN%*%t(zMN)) + rMN - invvMN <- MASS::ginv(vMN) - hMN <- MASS::ginv(t(xMN)%*%invvMN%*%xMN) - qMN <- invvMN - invvMN%*%xMN%*%(hMN)%*%t(xMN)%*%invvMN - eblueVn<-mfl@beta - eblupVn<-gammaMN%*%t(zMN)%*%invvMN%*%(yVn-xMN%*%eblueVn) ## equivalent de ranef(mfl) - rownames(eblupVn) <- colnames(zMN) - - - ## Calculs of matrices and vectors used in graph diagnosics --------------------------------------------- - - ## Marginal and individual predictions, residuals and variances - marpredVn <- xMN%*%eblueVn - marresVn <- yVn - marpredVn - marvarMN <- vMN - xMN%*%hMN%*%t(xMN) - - condpredVn <- marpredVn + zMN%*%eblupVn - condresVn <- yVn - condpredVn - condvarMN <- rMN%*%qMN%*%rMN - - - ## Analysis of marginal and conditional residuals - stmarresVn <-stcondresVn <- rep(0,nobsN) - lesverVn <- rep(0, nunitN) - names(lesverVn )<- idunitVc - - for(i in 1:nunitN){ - - idxiVn <- which(df[, unitC] == idunitVc[i]) ## position des observations du sujet i - miN <- length(idxiVn) - - ## standardization of marginal residual - stmarresVn[idxiVn] <- as.vector(solve(sqrtmF(marvarMN[idxiVn,idxiVn]))%*%marresVn[idxiVn]) - - ##Standardized Lessafre and Verbeke's measure - auxMN <- diag(1, ncol = miN, nrow =miN)- stmarresVn[idxiVn]%*%t(stmarresVn[idxiVn]) - lesverVn[i] <- sum(diag(auxMN%*%t(auxMN))) - - ## standardization of conditional residual - stcondresVn[idxiVn] <- as.vector(solve(sqrtmF(condvarMN[idxiVn,idxiVn]))%*%condresVn[idxiVn]) - } - lesverVn <- lesverVn/sum(lesverVn) - - - ## Least confounded residuals (à valider !) - ## résultats différents des valeurs estimées via le script de Singer, car - ## non unicité des vecteurs propres de la décomposition spectrale ? - if(least.confounded){ - sqrMN <- sqrtmF(rMN) - specDec <- eigen((sqrMN%*%qMN%*%sqrMN), symmetric = TRUE, only.values = FALSE) - - cMN <- t(sqrt(solve(diag((specDec$values[1:(nobsN -pN)])))) %*% t(specDec$vectors[1:(nobsN -pN),1:(nobsN -pN)]) %*% - solve(sqrtmF(rMN[1:(nobsN -pN),1:(nobsN -pN)])) ) - - lccondresVn <- (cMN%*%condresVn[1:(nobsN -pN)])/sqrt(diag(cMN%*%condvarMN[1:(nobsN -pN),1:(nobsN -pN)]%*%t(cMN))) - } - - - ## EBLUP analysis (Mahalanobis' distance) - varbMN <- gammaMN%*%t(zMN)%*%qMN%*%zMN%*%gammaMN - mdistVn <- rep(0, nunitN) - for(i in 1:nunitN){ - mdistVn[i] <- eblupVn[i]^2/varbMN[i, i] - } - mdistVn <- mdistVn/sum(mdistVn) - - - ## Combine data and results in 2 data frames for easy plotting with ggplot2 ----------------------------- - - ## long data frame (all observations) - - df <- data.frame(df, - mar.pred = marpredVn, - mar.res = marresVn, - st.mar.res = stmarresVn, - cond.pred = condpredVn, - cond.res = condresVn, - st.cond.res = stcondresVn) - - if(!("fixfact" %in% colnames(df))) - df$fixfact <- rep(1, nrow(df)) - df$numTime <- as.numeric(levels(as.factor(df$time))) - - - ## short data frame (1 row per unit) - - unitDf <- data.frame(unit = idunitVc, - eblup = eblupVn, - lvm = lesverVn, - mal = mdistVn) - - unitDf$fixfact <- sapply(1:nrow(unitDf), - function(i){ - unique(df[which(df[, unitC] == unitDf$unit[i]), - fixfactC]) - }) - - unitDf$se <- rep(NA, nrow(unitDf)) - re <- ranef(mfl, condVar =TRUE) - for(i in 1:nrow(unitDf)) - unitDf$se[i] <- sqrt(attr(re[[unitC]], "postVar")[1,1,i]) - unitDf$upper <- unitDf$eblup+unitDf$se*1.96 - unitDf$lower <- unitDf$eblup-unitDf$se*1.96 - - - ## Outliers "annotations" - df$marres.out <- rep("", nrow(df)) - df$marres.out[abs(df$st.mar.res)>hlimitN] <- - paste(df[abs(df$st.mar.res)>hlimitN, unitC], - df[abs(df$st.mar.res)>hlimitN, timeC], - sep = ".") - df$marres.out <- paste(" ", df$marres.out, sep ="") - - df$condres.out <- rep("", nrow(df)) - df$condres.out[abs(df$st.cond.res)>hlimitN] <- - paste(df[abs(df$st.cond.res)>hlimitN, unitC], - df[abs(df$st.cond.res)>hlimitN, timeC], - sep = ".") - df$condres.out <- paste(" ", df$condres.out, sep ="") - - ## Diagnostic Plots ------------------------------------------------------------------------------------- - - ## Linearity of effect and outlying observations - p1 <- ggplot(data = df, - aes(x=mar.pred, - y=st.mar.res, - colour=fixfact)) + - geom_point(size =2) + - geom_hline(yintercept = 0, col = "grey")+ - geom_smooth(aes(x=mar.pred, y=st.mar.res), data = df, se = FALSE, col = "blue", method = "loess")+ - ggtitle("Linearity of effects/outlying obervations")+ - xlab("Marginal predictions")+ - ylab("Standardized marginal residuals")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+ - geom_hline(yintercept = c(-1,1)*hlimitN, linetype = "dashed")+ - geom_text(aes(label = marres.out, col = fixfact), hjust=0, vjust=0) - - # p1hist <- histF(df$st.mar.res, sd_x =1)+ - # xlab("Standardized marginal residuals") - - - ## Presence of outlying observations and homoscedacity of residuals - p2 <- ggplot(data = df, - aes(x=cond.pred, - y=st.cond.res, - colour=fixfact)) + - geom_point(size =2) + - geom_hline(yintercept = 0, col = "grey")+ - geom_smooth(aes(x=cond.pred, y=st.cond.res), data = df, se = FALSE, col = "blue", method = "loess")+ - ggtitle("Homosedasticity of conditional residuals/outlying observations")+ - xlab("Individual predictions")+ - ylab("Standardized conditional residuals")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+ - geom_hline(yintercept = c(-1,1)*hlimitN, linetype = "dashed")+ - geom_text(aes(label = condres.out, col = fixfact), hjust=0, vjust=0) - - # p2hist <- histF(df$st.cond.res, sd_x =1)+ - # xlab("Standardized conditional residuals") - - - ## Normality of residuals - if(least.confounded){ - p3 <- qqplotF(x = lccondresVn, - distribution = "norm", - line.estimate = NULL, - conf = 0.95)+ - xlab("Standard normal quantiles")+ - ylab("Least confounded conditional residual quantiles")+ - ggtitle("Normality of conditional error (least confounded)")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold")) - - # p3hist <- histF(lccondresVn, sd_x =1)+ - # xlab("Least confounded conditional residuals") - }else{ - p3 <-qqplotF(x = df$st.cond.res, - distribution = "norm", - line.estimate = NULL, - conf = 0.95)+ - xlab("Standard normal quantiles")+ - ylab("Standardized conditional residual quantiles")+ - ggtitle("Normality of conditional error")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold")) - # p3hist <-p2hist - } - - ## Within-units covariance structure - p4 <- ggplot(data = unitDf, - aes(x=unit, - y=lvm, - colour=fixfact)) + - geom_point(size =2) + - theme(legend.position="none")+ - xlab(unitC)+ - ylab("Standardized Lesaffre-Verbeke measure")+ - geom_hline(yintercept = 2*mean(unitDf$lvm), linetype = "dashed")+ - geom_text(aes(label = unit, col = fixfact), - data = unitDf[unitDf$lvm>2*mean(unitDf$lvm), ], - hjust=0, vjust=0)+ - ggtitle("Within-units covariance matrice")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold")) - - ## EBLUP modif lmerTest v3 mais plante si pas d'effet aléatoire - #pvl <- ranova(model=mfl, reduce.terms = TRUE) - #pvrnd <- pvl[[6]][2] - #ggtitle(paste("Random effect on intercept (LRT p-value = ",round(pvrnd,digits = 5), ")", sep = ""))+ - - p5 <- - ggplot(aes(x = eblup, y = unit, col = fixfact), data = unitDf)+ - geom_point(size =3)+ - geom_segment(aes(xend = lower, yend = unit)) + - geom_segment(aes(xend = upper, yend = unit))+ - ggtitle("Random effect on intercept")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+ - geom_vline(xintercept = 0, linetype = "dashed")+ - ylab(unitC) - - ## p6 - p6 <-qqplotF(x = unitDf$mal, - distribution = "chisq", - df= 1, - line.estimate = NULL, - conf = 0.95)+ - xlab("Chi-squared quantiles")+ - ylab("Standadized Mahalanobis distance")+ - ggtitle("Normality of random effect")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold")) - - ## Outlying subjects - p7 <- - ggplot(aes(y = mal, x= unit, col = fixfact), data = unitDf)+ - geom_point(size =3)+ - ylab("Standardized Mahalanobis distance")+ - geom_vline(xintercept = 0, linetype = "dashed")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))+ - geom_hline(yintercept = 2*mean(unitDf$mal), linetype = "dashed")+ - geom_text(aes(label = unit, col = fixfact), - data = unitDf[unitDf$mal>2*mean(unitDf$mal), ], - hjust=1, vjust=0)+ - ggtitle("Outlying subjects")+ - xlab(unitC) - - ## "Data" and "modeling" Plots -------------------------------------------------------------------------- - - ## Individual time-course - rawPlot <- ggplot(data = df, - aes(x=time, y=vd, colour=fixfact, group=subject)) + - geom_line() + ggtitle("Individual time-courses ")+ - theme(legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold")) - - ## Post-hoc estimates (modification due to lmerTest v3) - ddlsm1 <- data.frame(difflsmeans(mfl,test.effs=NULL)) - colnames(ddlsm1)[9] <- "pvalue" - # ddlsm1$name <- sapply(rownames(ddlsm1), - # function(nam){ - # strsplit(nam, split = " ", fixed =TRUE)[[1]][1] - # }) - # ddlsm1$detail <- sapply(rownames(ddlsm1), - # function(nam){ - # paste(strsplit(nam, split = " ", fixed =TRUE)[[1]][-1], - # collapse= "") - # }) - # - # colnames(ddlsm1)<- make.names(colnames(ddlsm1)) - ddlsm1$Significance <- rep("NS", nrow(ddlsm1)) - - ## modif JF pour tenir compte du seuil de pvalues defini par le user - ddlsm1$Significance[which(ddlsm1$pvalue <pvalCutof & ddlsm1$pvalue >=0.01)] <- "p-value < threshold" - ddlsm1$Significance[which(ddlsm1$pvalue <0.01 & ddlsm1$pvalue >=0.005)] <- "p-value < 0.01" - ddlsm1$Significance[which(ddlsm1$pvalue <0.005)] <- "p-value < 0.005" - - phPlot <- ggplot(ddlsm1, aes(x = levels, y = Estimate))+ - facet_grid(facets = ~term, ddlsm1,scales = "free", space = "free")+ - geom_bar( aes(fill = Significance), stat="identity")+ - theme(axis.text.x = element_text(angle = 90, hjust = 1))+ - scale_fill_manual( - values = c("NS" = "grey", - "p-value < threshold" = "yellow", - "p-value < 0.01" = "orange", - "p-value < 0.005" = "red"))+ - geom_errorbar(aes(ymin = lower, ymax =upper ), width=0.25)+ - ggtitle("Post-hoc estimates")+xlab("")+ - theme(plot.title = element_text(size = rel(1.2), face = "bold")) - - ## Final plotting - - blank<-rectGrob(gp=gpar(col="white")) - rectspacer<-rectGrob(height = unit(0.1, "npc"), gp=gpar(col="grey")) - - grid.arrange(blank, - rawPlot, blank, phPlot, - rectspacer, - p1,blank, p2, blank, p3, blank, p4, - blank, - p5,blank, p6, blank, p7, blank, - blank,blank, - top = textGrob(title,gp=gpar(fontsize=40,font=4)), - layout_matrix = matrix(c(rep(1,7), - 2, 3, rep(4,3), 20,21, - rep(5,7), - 6:12, - rep(13,7), - 14:18, rep(19,2)), - ncol=7, nrow=6, byrow=TRUE), - heights= c(0.1/3, 0.3, 0.1/3, 0.3, 0.1/3, 0.3 ), - widths = c(0.22, 0.04, 0.22,0.04 , 0.22, 0.04, 0.22)) - - invisible(NULL) - -} - -#diagmflF(mfl, title = "", outlier.limit = 3, least.confounded =TRUE) - -plot.res.Lmixed <- function(mfl, df, title = "", pvalCutof = 0.05) { - - nameVar <- colnames(df)[4] - fflab <- colnames(df)[1] - ## Individual time-course - rawPlot <- - ggplot(data = df, aes(x=df[[2]], y=df[[4]], colour=df[[1]], group=df[[3]])) + - geom_line() + ggtitle("Individual time-courses (raw data)")+ - ylab(nameVar) + - xlab(label = colnames(df)[2])+ - theme(legend.title = element_blank() , legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold")) - - ## Boxplot of fixed factor - bPlot <- - ggplot(data = df, aes(y=df[[4]], x=df[[1]], color=df[[1]]))+ - geom_boxplot(outlier.colour="red", outlier.shape=8, outlier.size=4)+ - ggtitle(paste("Boxplot by ",fflab,sep=""))+xlab("")+ylab("")+ - theme(legend.title = element_blank(), plot.title = element_text(size = rel(1.2), face = "bold")) - - ## Post-hoc estimates - - ddlsm1 <- mfl - ddlsm1$name <- rownames(ddlsm1) - # ddlsm1$name <- sapply(rownames(ddlsm1), - # function(nam){ - # strsplit(nam, split = " ", fixed =TRUE)[[1]][1] - # }) - # ddlsm1$detail <- sapply(rownames(ddlsm1), - # function(nam){ - # paste(strsplit(nam, split = " ", fixed =TRUE)[[1]][-1], - # collapse= "") - # }) - # - #colnames(ddlsm1)<- make.names(colnames(ddlsm1)) - ddlsm1$Significance <- rep("NS", nrow(ddlsm1)) - - ## modif JF pour tenir compte du seuil de pvalues defini par le user - options("scipen"=100, "digits"=5) - pvalCutof <- as.numeric(pvalCutof) - bs=0.05; bm=0.01; bi=0.005 - if (pvalCutof >bm) {bs <- pvalCutof} else - if (pvalCutof <bm & pvalCutof >bi) {bm <- pvalCutof; bs <- pvalCutof} else - if (pvalCutof < bi) {bi <- pvalCutof; bm <- pvalCutof; bs <- pvalCutof} - lbs <- paste("p-value < ",bs, sep="") - lbm <- paste("p-value < ",bm, sep="") - lbi <- paste("p-value < ",bi, sep="") - - cols <- paste("p-value < ",bs, sep="") - colm <- paste("p-value < ",bm, sep="") - coli <- paste("p-value < ",bi, sep="") - valcol <- c("grey","yellow","orange","red") - names (valcol) <- c("NS",lbs,lbm,lbi) - ddlsm1$Significance[which(ddlsm1$p.value<= bs)] <- lbs - ddlsm1$Significance[which(ddlsm1$p.value<bs & ddlsm1$p.value>= bm)] <- lbm - ddlsm1$Significance[which(ddlsm1$p.value<bi)] <- lbi - - phPlot <- - ggplot(ddlsm1, aes(x = levels, y = Estimate))+ - facet_grid(facets = ~term, ddlsm1,scales = "free", space = "free")+ - geom_bar( aes(fill = Significance), stat="identity")+ - theme(axis.text.x = element_text(angle = 90, hjust = 1))+ - scale_fill_manual( - # values = c("NS" = "grey", "p-value < threshold" = "yellow","p-value < 0.01" = "orange","p-value < 0.005" = "red"))+ - # values = c("NS" = 'grey', "pvalue < 0.05 "= 'yellow',"p-value < 0.01" = 'orange',"p-value < 0.005" = 'red'))+ - # values = c("NS = grey", "p-value < threshold = yellow","p-value < 0.01 = orange","p-value < 0.005 = red"))+ - values = valcol )+ - geom_errorbar(aes(ymin = Lower.CI, ymax =Upper.CI ), width=0.25)+ - # ggtitle(paste("Post-hoc estimates with p-value threshold = ",pvalCutof,sep=""))+xlab("")+ - ggtitle("Post-hoc estimates ")+xlab("")+ - theme(plot.title = element_text(size = rel(1.2), face = "bold")) - - ## Final plotting - grid.arrange(arrangeGrob(rawPlot,bPlot,ncol=2), - phPlot,nrow=2, - top = textGrob(title,gp=gpar(fontsize=32,font=4)) - ) - -} \ No newline at end of file