| Previous changeset 0:1422de181204 (2018-10-10) |
|
Commit message:
Uploaded |
|
added:
mixmodel4repeated_measures/diagmfl.R mixmodel4repeated_measures/mixmodel.xml mixmodel4repeated_measures/mixmodel_script.R mixmodel4repeated_measures/mixmodel_wrapper.R mixmodel4repeated_measures/test-data/OneFactor_dataMatrix.tsv mixmodel4repeated_measures/test-data/OneFactor_sampleMetadata.tsv mixmodel4repeated_measures/test-data/OneFactor_variableMetadata.tsv mixmodel4repeated_measures/test-data/TwoFactor_dataMatrix.txt mixmodel4repeated_measures/test-data/TwoFactor_sampleMetadata.txt mixmodel4repeated_measures/test-data/TwoFactor_variableMetadata.txt mixmodel4repeated_measures/test-data/mixmodel_TwoFactor_variableMetadata.txt |
|
removed:
diagmfl.R mixmodel.xml mixmodel_script.R mixmodel_wrapper.R test-data/demo1_Samples.txt test-data/demo1_matrix.txt test-data/demo1_variables.txt test-data/mixmodel_demo1_variables.txt test-data/test_raw_variable_out.txt |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 diagmfl.R --- a/diagmfl.R Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,568 +0,0 @@\n-###############################################################################################\n-## Diagnostics graphics pour les mod\xe8les lin\xe9aires mixtes de type "mfl" ##\n-############################################################################################### ##\n-## Input: ##\n-## mfl, un mod\xe8le lin\xe9aire mixte g\xe9n\xe9r\xe9 par le module lmixedm de JF Martin (fonction ##\n-## lmRepeated2FF), i.e. : un mod\xe8le lin\xe9aire mixte de type lmer cr\xe9\xe9 par la formule ##\n-## mfl <- lmer( vd ~ time + fixfact + time:fixfact + (1| subject), ids) ##\n-## => noms de colonnes importants ##\n-## => 1 seul effet al\xe9atoire sur l\'intercept ##\n-###############################################################################################\n-## Output : ##\n-## Les graphics, les calculs associ\xe9s et les notations* utilis\xe9es dans le script suivent ##\n-## l\'article de Singer et al (2016) Graphical Tools for detedcting departures from linear ##\n-## mixed model assumptions and some remedial measures, International Statistical Review ##\n-## (doi:10.1111/insr.12178) ##\n-## * ajout d\'une ou 2 lettres de typage de la variables ##\n-##\n-## Script adapt\xe9 de http://www.ime.unicamp.br/~cnaber/residdiag_nlme_v22.R pour fonction- ##\n-## ner avec un mod\xe8le lmer (et non lme), des sujets avec des identifiants non num\xe9riques, ##\n-## et des observations non ordonn\xe9es sujet par sujet (dernier point \xe0 v\xe9rifier.) ##\n-## ############################################################################################\n-## Remarques sur les calculs num\xe9riques ##\n-## - l\'inverse d\'une matrice est calcul\xe9e \xe0 partir de la fonction ginv du package MASS ##\n-## (Moore- Penrose generalized inverse) au lieu de la fonction "solve" ##\n-## - la racine carr\xe9e des matrices sont calcul\xe9es gr\xe2ce \xe0 une d\xe9omposition SVD (car ##\n-## s\'applique normalement ici uniquement \xe0 des matrices sym\xe9triques positives). A ##\n-## remplacer par la fonction sqrtm{expm} si des erreurs apparaissent ??? ##\n-###############################################################################################\n-\n-\n-library(ggplot2)\n-library(gridExtra)\n-library(grid)\n-\n-\n-##-------------------------------------------------------------------------------------------------##\n-## Helpers\n-##-------------------------------------------------------------------------------------------------##\n-\n-## square root of a matrix\n-## http://www.cs.toronto.edu/~jepson/csc420/notes/introSVD.pdf (page 6)\n-## (for matMN a n x n matrix that symetric and non-negative definite)\n-\n-sqrtmF<-function(matMN){\n- matMN <- as.matrix(matMN)\n- # ## check that matMN is symetric\n- # if(!all(t(matMN==matMN)))\n- # stop("matMN must be symetric.")\n- svd_dec <- svd(matMN)\n- invisible(svd_dec$u%*%sqrt(diag(svd_dec$d))%*%t(svd_dec$v))\n-}\n-\n-\n-## qqplotF\n-## adapted from https://gist.github.com/rentrop/d39a8406ad8af2a1066c\n-\n-\n-qqplotF <- function(x, \n- distribution = "norm", ..., \n- line.estimate = NULL, \n- conf = 0.95,\n- labels = names(x)){\n- q.function <- eval(parse(text = paste0("q", distribution)))\n- d.function <- eval(parse(text = paste0("d", distribution)))\n- x <- na.omit(x)\n- ord <- order(x)\n- n <- length(x)\n- P <- ppoints(length(x))\n- daf <- data.frame(ord.x = x[ord], z = q.function('..b' ncol=7, nrow=6, byrow=TRUE),\n- heights= c(0.1/3, 0.3, 0.1/3, 0.3, 0.1/3, 0.3 ),\n- widths = c(0.22, 0.04, 0.22,0.04 , 0.22, 0.04, 0.22))\n- \n- invisible(NULL)\n-\n-}\n-\n-#diagmflF(mfl, title = "", outlier.limit = 3, least.confounded =TRUE)\n-\n-plot.res.Lmixed <- function(mfl, df, title = "", pvalCutof = 0.05) {\n-\n- nameVar <- colnames(df)[4]\n- fflab <- colnames(df)[1]\n- ## Individual time-course\n- rawPlot <- \n- ggplot(data = df, aes(x=df[[2]], y=df[[4]], colour=df[[1]], group=df[[3]])) +\n- geom_line() + ggtitle("Individual time-courses (raw data)")+ \n- ylab(nameVar) +\n- xlab(label = colnames(df)[2])+\n- theme(legend.title = element_blank() , legend.position="none", plot.title = element_text(size = rel(1.2), face = "bold"))\n- \n- ## Boxplot of fixed factor\n- bPlot <- \n- ggplot(data = df, aes(y=df[[4]], x=df[[1]], color=df[[1]]))+\n- geom_boxplot(outlier.colour="red", outlier.shape=8, outlier.size=4)+\n- ggtitle(paste("Boxplot by ",fflab,sep=""))+xlab("")+ylab("")+\n- theme(legend.title = element_blank(), plot.title = element_text(size = rel(1.2), face = "bold"))\n- \n- ## Post-hoc estimates\n- \n- ddlsm1 <- mfl\n- ddlsm1$name <- rownames(ddlsm1)\n- # ddlsm1$name <- sapply(rownames(ddlsm1),\n- # function(nam){\n- # strsplit(nam, split = " ", fixed =TRUE)[[1]][1]\n- # })\n- # ddlsm1$detail <- sapply(rownames(ddlsm1),\n- # function(nam){\n- # paste(strsplit(nam, split = " ", fixed =TRUE)[[1]][-1],\n- # collapse= "")\n- # })\n- # \n- #colnames(ddlsm1)<- make.names(colnames(ddlsm1))\n- ddlsm1$Significance <- rep("NS", nrow(ddlsm1))\n- \n- ## modif JF pour tenir compte du seuil de pvalues defini par le user \n- options("scipen"=100, "digits"=5)\n- pvalCutof <- as.numeric(pvalCutof)\n- bs=0.05; bm=0.01; bi=0.005\n- if (pvalCutof >bm) {bs <- pvalCutof} else\n- if (pvalCutof <bm & pvalCutof >bi) {bm <- pvalCutof; bs <- pvalCutof} else\n- if (pvalCutof < bi) {bi <- pvalCutof; bm <- pvalCutof; bs <- pvalCutof}\n- lbs <- paste("p-value < ",bs, sep="")\n- lbm <- paste("p-value < ",bm, sep="")\n- lbi <- paste("p-value < ",bi, sep="")\n- \n- cols <- paste("p-value < ",bs, sep="")\n- colm <- paste("p-value < ",bm, sep="")\n- coli <- paste("p-value < ",bi, sep="")\n- valcol <- c("grey","yellow","orange","red")\n- names (valcol) <- c("NS",lbs,lbm,lbi)\n- ddlsm1$Significance[which(ddlsm1$p.value<= bs)] <- lbs \n- ddlsm1$Significance[which(ddlsm1$p.value<bs & ddlsm1$p.value>= bm)] <- lbm\n- ddlsm1$Significance[which(ddlsm1$p.value<bi)] <- lbi\n-\n- phPlot <- \n- ggplot(ddlsm1, aes(x = levels, y = Estimate))+\n- facet_grid(facets = ~term, ddlsm1,scales = "free", space = "free")+\n- geom_bar( aes(fill = Significance), stat="identity")+\n- theme(axis.text.x = element_text(angle = 90, hjust = 1))+\n- scale_fill_manual(\n- # values = c("NS" = "grey", "p-value < threshold" = "yellow","p-value < 0.01" = "orange","p-value < 0.005" = "red"))+\n- # values = c("NS" = \'grey\', "pvalue < 0.05 "= \'yellow\',"p-value < 0.01" = \'orange\',"p-value < 0.005" = \'red\'))+\n- # values = c("NS = grey", "p-value < threshold = yellow","p-value < 0.01 = orange","p-value < 0.005 = red"))+\n- values = valcol )+ \n- geom_errorbar(aes(ymin = Lower.CI, ymax =Upper.CI ), width=0.25)+\n- # ggtitle(paste("Post-hoc estimates with p-value threshold = ",pvalCutof,sep=""))+xlab("")+\n- ggtitle("Post-hoc estimates ")+xlab("")+\n- theme(plot.title = element_text(size = rel(1.2), face = "bold"))\n- \n- ## Final plotting\n- grid.arrange(arrangeGrob(rawPlot,bPlot,ncol=2),\n- phPlot,nrow=2,\n- top = textGrob(title,gp=gpar(fontsize=32,font=4))\n- )\n- \n-} \n\\ No newline at end of file\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel.xml --- a/mixmodel.xml Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,211 +0,0 @@\n-<tool id="mixmodel" name="mixmodel" version="1.0.0">\n- <description>ANOVA for repeated measures statistics</description>\n-\n- <requirements>\n- <requirement type="package" version="1.1_4">r-batch</requirement>\n- <requirement type="package" version="1.1_13">r-lme4</requirement>\n- <requirement type="package" version="3.0_1">r-lmertest</requirement>\n- <requirement type="package" version="2.34.0">bioconductor-multtest</requirement>\n- <requirement type="package" version="2.3">r-gridextra</requirement>\n- </requirements>\n-\n- <stdio>\n- <exit_code range="1:" level="fatal" />\n- </stdio>\n-\n- <command><![CDATA[\n- Rscript $__tool_directory__/mixmodel_wrapper.R\n-\n- dataMatrix_in "$dataMatrix_in"\n- sampleMetadata_in "$sampleMetadata_in"\n- variableMetadata_in "$variableMetadata_in"\n-\n- fixfact "$fixfact"\n- time "$time"\n- subject "$subject"\n- adjC "$adjC"\n- trf "$trf"\n- thrN "$thrN"\n- diaR "$diaR"\n-\n- variableMetadata_out "$variableMetadata_out"\n- out_graph_pdf "$out_graph_pdf"\n-\t\tout_estim_pdf "$out_estim_pdf"\n- information "$information"\n-\n- ]]></command>\n-\n- <inputs>\n- <param name="dataMatrix_in" label="Data matrix file" type="data" format="tabular" help="variable x sample, decimal: \'.\', missing: NA, mode: numerical, sep: tabular" />\n- <param name="sampleMetadata_in" label="Sample metadata file" type="data" format="tabular" help="sample x metadata, decimal: \'.\', missing: NA, mode: character and numerical, sep: tabular" />\n- <param name="variableMetadata_in" label="Variable metadata file" type="data" format="tabular" help="variable x metadata, decimal: \'.\', missing: NA, mode: character and numerical, sep: tabular" />\n- <param name="fixfact" label="Fixed Factor of interest" type="text" help="Name of the column of the sample metadata table corresponding to the fixed factor"/>\n- <param name="time" label="Repeated factor (time)" type="text" help="Name of the column of the sample metadata table corresponding to the repeated factor"/>\n- <param name="subject" label="Subject factor" type="text" help="Name of the column of the sample metadata table corresponding to the subject factor"/>\n- <param name="adjC" label="Method for multiple testing correction" type="select" help="">\n- <option value="fdr">fdr</option>\n- <option value="BH">BH</option>\n- <option value="bonferroni">bonferroni</option>\n- <option value="BY">BY</option>\n- <option value="hochberg">hochberg</option>\n- <option value="holm">holm</option>\n- <option value="hommel">hommel</option>\n- <option value="none">none</option>\n- </param>\n- <param name="trf" label="Log transform of raw data" type="select" help="Transformation of raw data">\n- <option value="none">none</option>\n- <option value="log10">log10</option>\n- <option value="log2">log2</option>\n- </param>\n- <param name="thrN" type="float" value="0.05" label="(Corrected) p-value significance threshold" help="Must be between 0 and 1"/>\n- <param name="diaR" label="Perform diagnostic of the residuals" type="select"\n- help=" Used to assess the quality of models considering distribution of residuals ">\n- <option value="yes">yes</option>\n- <option value="no"></option>\n- </param>\n-\n- </inputs>\n-\n- <outputs>\n- <data name="variableMetadata_out" label="${tool.name}_${variableMetadata_in.name}" format="tabular"/>\n- <data name="information" label="${tool.name}_information.txt" format="txt"/>\n- <data name="out_graph_pdf" label="${tool.name}_diagResiduals" format="pdf"/>\n- <data name="out_estim_pdf" label="${tool.name}_Estimates" format="pdf"/>\n-\n- </outputs>\n-\n- <tests>\n- <test>\n- '..b' the table must not contain metadata apart from row and column names; the row and column names must be identical to the rownames of the sample and variable metadata, respectively (see below)\n-|\n-\n-Sample metadata file\n-| sample x metadata **sampleMetadata** tabular separated file of the numeric and/or character sample metadata, with . as decimal and NA for missing values\n-|\n-\n-Variable metadata file\n-| variable x metadata **variableMetadata** tabular separated file of the numeric and/or character variable metadata, with . as decimal and NA for missing values\n-|\n-\n-\n-Treatment\n-| Name of the fixed factor in the sample metadata file\n-|\n-\n-Time\n-| Name of the repeated (time) factor in the sample metadata\tfile\n-|\n-\n-Subject\n-| Name of the subject (on which the repeated measurement id done) in the sample metadata file\n-|\n-\n-Method for multiple testing correction\n-| The 7 methods implemented in the \'p.adjust\' R function are available and documented as follows:\n-| "The adjustment methods include the Bonferroni correction ("bonferroni") in which the p-values are multiplied by the number of comparisons. Less conservative corrections are also included by Holm (1979) ("holm"), Hochberg (1988) ("hochberg"), Hommel (1988) ("hommel"), Benjamini and Hochberg (1995) ("BH" or its alias "fdr"), and Benjamini and Yekutieli (2001) ("BY"), respectively. A pass-through option ("none") is also included. The set of methods are contained in the p.adjust.methods vector for the benefit of methods that need to have the method as an option and pass it on to p.adjust. The first four methods are designed to give strong control of the family-wise error rate. There seems no reason to use the unmodified Bonferroni correction because it is dominated by Holm\'s method, which is also valid under arbitrary assumptions. Hochberg\'s and Hommel\'s methods are valid when the hypothesis tests are independent or when they are non-negatively associated (Sarkar, 1998; Sarkar and Chang, 1997). Hommel\'s method is more powerful than Hochberg\'s, but the difference is usually small and the Hochberg p-values are faster to compute. The "BH" (aka "fdr") and "BY" method of Benjamini, Hochberg, and Yekutieli control the false discovery rate, the expected proportion of false discoveries amongst the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so these methods are more powerful than the others."\n-\n-\n-(Corrected) p-value significance threshold\n-|\n-|\n-\n-------------\n-Output files\n-------------\n-\n-variableMetadata_out.tabular\n-| **variableMetadata** file identical to the file given as argument plus\n-| pvalue of Shapiro normality test of the residuals\n-| pvalues of the main effects and interaction\n-| PostHoc test with difference between levels and pvalues of these difference\n-|\n-\n-mixedmodel_diagResiduals\n-| if Perform diagnostic of the residuals" is set to yes(default) a pdf file is created with a graphical\n-| representation of differences among levels of factors with a color code for significance and an error bar\n-| Then a serie of graphics are output in order to assess the distribution of residuals to check the adjustment.\n-\n-information.txt\n-| File with all messages and warnings generated during the computation\n-| The list of variables with name and a tag if it is significant for at least fixed or repeated factor.\n-\n-\n- ]]></help>\n-\n- <citations>\n- <citation type="doi">10.18637/jss.v082.i13.</citation>\n- <citation type="bibtex">@ARTICLE{fisher,\n- author = {Benjamini Y. and Hochberg Y.,\n- title = {Controlling the false discovery rate: a practical and powerful approach for multiple testing. Journal of the Royal Statistical Society},\n- journal = {Series B (Methodological)},\n- year = {1995},\n- volume = {57},\n- pages = {289-300}\n- }</citation>\n- <citation type="doi">10.1093/bioinformatics/btu813</citation>\n- </citations>\n-\n-</tool>\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/diagmfl.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/diagmfl.R Mon May 16 12:31:58 2022 +0000 |
| [ |
| b'@@ -0,0 +1,918 @@\n+#\' Calcul des grandeurs "diagnostiques"\r\n+#\'\r\n+#\' Script adapte de http://www.ime.unicamp.br/~cnaber/residdiag_nlme_v22.R pour fonctionner\r\n+#\' avec un modele lmer (et non lme), des sujets avec des identifiants non numeriques,\r\n+#\' et des observations non ordonnees sujet par sujet (dernier point a verifier.)\r\n+#\'\r\n+#\' @detail Les graphiques, les calculs associ\xc3\xa9s et les notations utilisees dans le script suivent\r\n+#\' l\'article de Singer et al (2016) Graphical Tools for detedcting departures from linear\r\n+#\' mixed model assumptions and some remedial measures, International Statistical Review\r\n+#\' (doi:10.1111/insr.12178)\r\n+#\'\r\n+#\' @param mfl A linear mixed model fitted via lmer or a data frame containing data\r\n+#\' @return A list\r\n+#\' @author Natacha Lenuzza\r\n+#\' @examples\r\n+#\' print("hello !")\r\n+#\'\r\n+#\' @export lmer.computeDiag\r\n+\r\n+\r\n+\r\n+lmer.computeDiag <- function(mfl) {\r\n+\r\n+ ## Check arguments ---------------------------------------------------------\r\n+\r\n+ if (length(mfl@flist) > 1)\r\n+ stop("Several \'grouping level\' for random effect not implemented yet.")\r\n+\r\n+\r\n+ ## extracting information from mfl models -------------------------------------------------------------\r\n+ # data\r\n+ df <- mfl@frame\r\n+ responseC <- names(df)[1]\r\n+ unitC <- names(mfl@flist)[1]\r\n+ # observations\r\n+ yVn <- df[, responseC]\r\n+ nobsN <- length(yVn)\r\n+ # units\r\n+ idunitVc <- levels(mfl@flist[[1]])\r\n+ nunitN <- length(unique(idunitVc))\r\n+ #X\r\n+ xMN <- mfl@pp$X\r\n+ pN <- ncol(xMN)\r\n+ #Z\r\n+ zMN <- t(as.matrix(mfl@pp$Zt))\r\n+ # Estimated covariance matrix of random effects (Gam)\r\n+ aux <- VarCorr(mfl)[[1]] ## assuming only one level of grouping\r\n+ aux2 <- attr(aux, "stddev")\r\n+ gMN <- attr(aux, "correlation") * (aux2 %*% t(aux2))\r\n+ gammaMN <- as.matrix(kronecker(diag(nunitN), gMN))\r\n+ q <- dim(gMN)[1]\r\n+ # Estimated covariance matrix of conditonal error (homoskedastic conditional independance model)\r\n+ sigsqN <- attr(VarCorr(mfl), "sc")^2\r\n+ rMN <- sigsqN * diag(nobsN)\r\n+ # Estimated covariance matrix of Y\r\n+ vMN <- (zMN %*% gammaMN %*% t(zMN)) + rMN\r\n+ invvMN <- MASS::ginv(vMN)\r\n+ # H and Q matrix\r\n+ hMN <- MASS::ginv(t(xMN) %*% invvMN %*% xMN)\r\n+ qMN <- invvMN - invvMN %*% xMN %*% (hMN) %*% t(xMN) %*% invvMN\r\n+ # eblue and eblup\r\n+ eblueVn <- mfl@beta\r\n+ eblupVn <- gammaMN %*% t(zMN) %*% invvMN %*% (yVn - xMN %*% eblueVn) ## equivalent de ranef(mfl)\r\n+ rownames(eblupVn) <- colnames(zMN)\r\n+ ## Calculs of matrices and vectors used in graph diagnosics ---------------------------------------------\r\n+ ## Marginal and individual predictions, residuals and variances\r\n+ marpredVn <- xMN %*% eblueVn\r\n+ marresVn <- yVn - marpredVn\r\n+ marvarMN <- vMN - xMN %*% hMN %*% t(xMN)\r\n+ condpredVn <- marpredVn + zMN %*% eblupVn\r\n+ condresVn <- yVn - condpredVn\r\n+ condvarMN <- rMN %*% qMN %*% rMN\r\n+ ## Analysis of marginal and conditional residuals\r\n+ stmarresVn <- stcondresVn <- rep(0, nobsN)\r\n+ lesverVn <- rep(0, nunitN)\r\n+ names(lesverVn) <- idunitVc\r\n+ for (i in 1:nunitN) {\r\n+ idxiVn <- which(df[, unitC] == idunitVc[i]) ## position des observations du sujet i\r\n+ miN <- length(idxiVn)\r\n+ ## standardization of marginal residual\r\n+ stmarresVn[idxiVn] <- as.vector(solve(sqrtmF(marvarMN[idxiVn, idxiVn])) %*% marresVn[idxiVn])\r\n+ ##Standardized Lessafre and Verbeke\'s measure\r\n+ auxMN <- diag(1, ncol = miN, nrow = miN) - stmarresVn[idxiVn] %*% t(stmarresVn[idxiVn])\r\n+ lesverVn[i] <- sum(diag(auxMN %*% t(auxMN)))\r\n+ ## standardization of conditional residual\r\n+ stcondresVn[idxiVn] <- as.vector(solve(sqrtmF(condvarMN[idxiVn, idxiVn])) %*% condresVn[idxiVn])\r\n+ }\r\n+ lesverVn <- lesverVn / sum(lesverVn)\r\n+ ## Least confounded conditional residuals\r\n+ ## EBLUP analysis (Mahalanobis\' distance)\r\n+ varbMN <- gammaMN %*% t(zMN) %*% qMN %*% zMN %*% gammaMN\r\n+ mdistVn <- rep(0, nunitN)\r\n+ qm <- q - 1\r\n+ for (j in 1:nunitN) {\r\n+ gbi <- varbMN[(q * j - qm):'..b'diff(range(x)) / (2 * IQR(x) / length(x) ^ (1 / 3))\r\n+ )\r\n+ }\r\n+\r\n+\r\n+ daf <- data.frame(x = x)\r\n+ ## graph\r\n+ return(ggplot(data = daf, aes(x)) +\r\n+ geom_histogram(aes(y = ..density..),\r\n+ col = "black", fill = "grey", binwidth = bw) +\r\n+ geom_density(size = 1.2,\r\n+ col = "blue",\r\n+ linetype = "blank",\r\n+ fill = rgb(0, 0, 1, 0.1)) +\r\n+ stat_function(fun = dnorm,\r\n+ args = list(mean = 0, sd = sd_x),\r\n+ col = "blue", size = 1.2) +\r\n+ theme(legend.position = "none") +\r\n+ xlab(""))\r\n+}\r\n+\r\n+\r\n+plot.res.Lmixed <- function(mfl, df, title = "", pvalCutof = 0.05) {\r\n+\r\n+ ## define subscript of the different columns depending if we have only time (ncol(df)=3) or not\r\n+ if (ncol(df) > 3) {\r\n+ varidx <- 4\r\n+ ffidx <- 1\r\n+ timidx <- 2\r\n+ individx <- 3\r\n+ } else {\r\n+ varidx <- 3\r\n+ ffidx <- 1\r\n+ timidx <- 1\r\n+ individx <- 2\r\n+ }\r\n+ nameVar <- colnames(df)[varidx]\r\n+ fflab <- colnames(df)[ffidx]\r\n+ ## Individual time-course\r\n+ rawPlot <-\r\n+ ggplot(data = df, aes(x = df[[timidx]], y = df[[varidx]], colour = df[[ffidx]], group = df[[individx]])) +\r\n+ geom_point() +\r\n+ geom_line() + ggtitle("Individual time-courses (raw data)") +\r\n+ ylab(nameVar) +\r\n+ xlab(label = colnames(df)[2]) +\r\n+ theme(legend.title = element_blank(), legend.position = "none", plot.title = element_text(size = rel(1.2), face = "bold"))\r\n+ ## Boxplot of fixed factor\r\n+ bPlot <-\r\n+ ggplot(data = df, aes(y = df[[varidx]], x = df[[ffidx]], color = df[[ffidx]])) +\r\n+ geom_boxplot(outlier.colour = "red", outlier.shape = 8, outlier.size = 4) +\r\n+ ggtitle(paste("Boxplot by ", fflab, sep = "")) + xlab("") + ylab("") +\r\n+ theme(legend.title = element_blank(), plot.title = element_text(size = rel(1.2), face = "bold"))\r\n+ ## Post-hoc estimates\r\n+ ddlsm1 <- mfl\r\n+ ddlsm1$name <- rownames(ddlsm1)\r\n+ ddlsm1$Significance <- rep("NS", nrow(ddlsm1))\r\n+ ## modif JF pour tenir compte du seuil de pvalues defini par le user\r\n+ options("scipen" = 100, "digits" = 5)\r\n+ pvalCutof <- as.numeric(pvalCutof)\r\n+ bs <- 0.05; bm <- 0.01; bi <- 0.005\r\n+ if (pvalCutof > bm) {\r\n+ bs <- pvalCutof\r\n+ } else\r\n+ if (pvalCutof < bm & pvalCutof > bi) {\r\n+ bm <- pvalCutof; bs <- pvalCutof\r\n+ } else\r\n+ if (pvalCutof < bi) {\r\n+ bi <- pvalCutof; bm <- pvalCutof; bs <- pvalCutof\r\n+ }\r\n+ lbs <- paste("p-value < ", bs, sep = "")\r\n+ lbm <- paste("p-value < ", bm, sep = "")\r\n+ lbi <- paste("p-value < ", bi, sep = "")\r\n+ cols <- paste("p-value < ", bs, sep = "")\r\n+ colm <- paste("p-value < ", bm, sep = "")\r\n+ coli <- paste("p-value < ", bi, sep = "")\r\n+ valcol <- c("grey", "yellow", "orange", "red")\r\n+ names(valcol) <- c("NS", lbs, lbm, lbi)\r\n+ ddlsm1$Significance[which(ddlsm1$p.value <= bs)] <- lbs\r\n+ ddlsm1$Significance[which(ddlsm1$p.value < bs & ddlsm1$p.value >= bm)] <- lbm\r\n+ ddlsm1$Significance[which(ddlsm1$p.value < bi)] <- lbi\r\n+ ddlsm1$levels <- rownames(ddlsm1)\r\n+ ddlsm1$term <- sapply(rownames(ddlsm1), function(namC) {\r\n+ strsplit(namC, split = " ", fixed = TRUE)[[1]][1]\r\n+ })\r\n+\r\n+ phPlot <-\r\n+ ggplot(ddlsm1, aes(x = levels, y = Estimate)) +\r\n+ facet_grid(facets = ~term, ddlsm1, scales = "free", space = "free") +\r\n+ geom_bar(aes(fill = Significance), stat = "identity") +\r\n+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +\r\n+ scale_fill_manual(\r\n+ values = valcol) +\r\n+ geom_errorbar(aes(ymin = Lower.CI, ymax = Upper.CI), width = 0.25) +\r\n+ ggtitle("Post-hoc estimates ") + xlab("") +\r\n+ theme(plot.title = element_text(size = rel(1.2), face = "bold"))\r\n+\r\n+ ## Final plotting\r\n+ grid.arrange(arrangeGrob(rawPlot, bPlot, ncol = 2),\r\n+ phPlot, nrow = 2,\r\n+ top = textGrob(title, gp = gpar(fontsize = 32, font = 4))\r\n+ )\r\n+\r\n+}\r\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/mixmodel.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/mixmodel.xml Mon May 16 12:31:58 2022 +0000 |
| [ |
| b'@@ -0,0 +1,264 @@\n+<tool id="mixmodel4repeated_measures" name="mixmodel" version="3.1.0" profile="19.01">\r\n+ <description>ANOVA for repeated measures statistics</description>\r\n+\r\n+ <requirements>\r\n+ <requirement type="package" version="1.1_27">r-lme4</requirement>\r\n+ <requirement type="package" version="3.1_3">r-lmertest</requirement>\r\n+ <requirement type="package" version="3.2_0">r-ggplot2</requirement>\r\n+ <requirement type="package" version="2.3">r-gridExtra</requirement>\r\n+ <requirement type="package" version="2.40.0">bioconductor-multtest</requirement>\r\n+ </requirements>\r\n+\r\n+ <command detect_errors="exit_code"><![CDATA[\r\n+ Rscript \'$__tool_directory__/mixmodel_wrapper.R\'\r\n+\r\n+ dataMatrix_in \'$dataMatrix_in\'\r\n+ sampleMetadata_in \'$sampleMetadata_in\'\r\n+ variableMetadata_in \'$variableMetadata_in\'\r\n+\r\n+ fixfact \'$fixfact\'\r\n+ time \'$time\'\r\n+ subject \'$subject\'\r\n+ adjC \'$adjC\'\r\n+ trf \'$trf\'\r\n+ thrN \'$thrN\'\r\n+ diaR \'$diaR\'\r\n+ dff \'$dff\'\r\n+ rounding \'${roundchoice.rounding}\'\r\n+ #if $roundchoice.rounding == \'yes\' :\r\n+ decplaces \'${roundchoice.decplaces}\'\r\n+ #end if\r\n+\r\n+ variableMetadata_out \'$variableMetadata_out\'\r\n+ out_graph_pdf \'$out_graph_pdf\'\r\n+ out_estim_pdf \'$out_estim_pdf\'\r\n+ information \'$information\'\r\n+\r\n+ ]]></command>\r\n+\r\n+ <inputs>\r\n+ <param name="dataMatrix_in" label="Data matrix file" type="data" format="tabular" help="variable x sample, decimal: \'.\', missing: NA, mode: numerical, sep: tabular" />\r\n+ <param name="sampleMetadata_in" label="Sample metadata file" type="data" format="tabular" help="sample x metadata, decimal: \'.\', missing: NA, mode: character and numerical, sep: tabular" />\r\n+ <param name="variableMetadata_in" label="Variable metadata file" type="data" format="tabular" help="variable x metadata, decimal: \'.\', missing: NA, mode: character and numerical, sep: tabular" />\r\n+ <param name="fixfact" label="Fixed Factor of interest" type="text" help="Name of sample metadata column corresponding to the fixed factor (use none if only time factor"/>\r\n+ <param name="time" label="Repeated factor (time)" type="text" help="Name of the column of the sample metadata table corresponding to the repeated factor"/>\r\n+ <param name="subject" label="Subject factor" type="text" help="Name of the column of the sample metadata table corresponding to the subject factor"/>\r\n+ <param name="dff" label="Degrees of freedom method for the F-tests" type="select" help="">\r\n+ <option value="Satt">Satterthwaite</option>\r\n+ <option value="KenR">Kenward-Roger</option>\r\n+ </param>\r\n+ <param name="adjC" label="Method for multiple testing correction" type="select" help="">\r\n+ <option value="fdr">fdr</option>\r\n+ <option value="BH">BH</option>\r\n+ <option value="bonferroni">bonferroni</option>\r\n+ <option value="BY">BY</option>\r\n+ <option value="hochberg">hochberg</option>\r\n+ <option value="holm">holm</option>\r\n+ <option value="hommel">hommel</option>\r\n+ <option value="none">none</option>\r\n+ </param>\r\n+ <param name="trf" label="Log transform of raw data" type="select" help="Transformation of raw data">\r\n+ <option value="none">none</option>\r\n+ <option value="log10">log10</option>\r\n+ <option value="log2">log2</option>\r\n+ </param>\r\n+ <param name="thrN" type="float" value="0.05" label="(Corrected) p-value significance threshold" help="Must be between 0 and 1"/>\r\n+ <param name="diaR" label="Perform diagnostic of the residuals" type="select" help=" Used to assess the quality of models considering distribution of residuals ">\r\n+ <option value="yes">yes</option>\r\n+ <option value="no">no</option>\r\n'..b'88) ("hommel"), Benjamini and Hochberg (1995) ("BH" or its alias "fdr"), and Benjamini and Yekutieli (2001) ("BY"), respectively. A pass-through option ("none") is also included. The set of methods are contained in the p.adjust.methods vector for the benefit of methods that need to have the method as an option and pass it on to p.adjust. The first four methods are designed to give strong control of the family-wise error rate. There seems no reason to use the unmodified Bonferroni correction because it is dominated by Holm\'s method, which is also valid under arbitrary assumptions. Hochberg\'s and Hommel\'s methods are valid when the hypothesis tests are independent or when they are non-negatively associated (Sarkar, 1998; Sarkar and Chang, 1997). Hommel\'s method is more powerful than Hochberg\'s, but the difference is usually small and the Hochberg p-values are faster to compute. The "BH" (aka "fdr") and "BY" method of Benjamini, Hochberg, and Yekutieli control the false discovery rate, the expected proportion of false discoveries amongst the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so these methods are more powerful than the others."\r\n+\r\n+\r\n+(Corrected) p-value significance threshold\r\n+|\r\n+|\r\n+\r\n+Rounding the result\'s numerical columns\r\n+| This parameter enables the choice of a given number of decimal places to be used to display results\' values in the output columns of the variableMetadata that are generated by the tool.\r\n+| If set to "Yes", you can indicate the number of decimal places wanted, and all the results\' values will be rounded according to this number.\r\n+|\r\n+\r\n+------------\r\n+Output files\r\n+------------\r\n+\r\n+variableMetadata_out.tabular\r\n+| **variableMetadata** file identical to the file given as argument plus\r\n+| pvalue of Shapiro normality test of the residuals\r\n+| pvalues of the main effects and interaction\r\n+| PostHoc test with difference between levels and pvalues of these difference\r\n+|\r\n+\r\n+mixedmodel_Estimates\r\n+| A pdf file is created with a graphical representation of differences among levels of factors with a color code for significance and an error bar.\r\n+| Plots are only provided for features with at least one significant factor. \r\n+\r\n+\r\n+mixedmodel_diagResiduals\r\n+| if Perform diagnostic of the residuals" is set to yes(default) a pdf file is created with a a serie of\r\n+| graphics outputed in order to assess the distribution of residuals to check the adjustment.\r\n+| Plots are only provided for features with at least one factor being significant based on the non-corrected p-values. \r\n+\r\n+information.txt\r\n+| File with all messages and warnings generated during the computation\r\n+| The list of variables with name and a tag if it is significant for at least fixed or repeated factor.\r\n+\r\n+\r\n+\r\n+---------------------------------------------------\r\n+\r\n+Changelog/News\r\n+--------------\r\n+\r\n+**Version 3.1.0 - April 2022**\r\n+\r\n+- NEW FEATURE: the tool now can be used without fixed factor by specify "none" in the concerned parameter field.\r\n+\r\n+- NEW FEATURE: addition of the possibility to choose the degrees of freedom computation method (previously using only the Satterthwaite method).\r\n+\r\n+- NEW FEATURE: addition of the possibility to choose the number of decimal places in the variable-metadata numerical columns\' output.\r\n+\r\n+\r\n+ ]]></help>\r\n+\r\n+ <citations>\r\n+ <citation type="doi">10.18637/jss.v082.i13.</citation>\r\n+ <citation type="bibtex">@ARTICLE{fisher,\r\n+ author = {Benjamini Y. and Hochberg Y.,\r\n+ title = {Controlling the false discovery rate: a practical and powerful approach for multiple testing. Journal of the Royal Statistical Society},\r\n+ journal = {Series B (Methodological)},\r\n+ year = {1995},\r\n+ volume = {57},\r\n+ pages = {289-300}\r\n+ }</citation>\r\n+ <citation type="doi">10.1093/bioinformatics/btu813</citation>\r\n+ </citations>\r\n+\r\n+</tool>\r\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/mixmodel_script.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/mixmodel_script.R Mon May 16 12:31:58 2022 +0000 |
| [ |
| b'@@ -0,0 +1,481 @@\n+####### R functions to perform linear mixed model for repeated measures\r\n+####### on a multi var dataset using 3 files as used in W4M\r\n+##############################################################################################################\r\n+lmRepeated2FF <- function(ids, ifixfact, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]],\r\n+ pvalCutof = 0.05, dffOption, visu = visu, tit = "", least.confounded = FALSE, outlier.limit = 3) {\r\n+ ### function to perform linear mixed model with 1 Fixed factor + Time + random factor subject\r\n+ ### based on lmerTest package providing functions giving the same results as SAS proc mixed\r\n+ options(scipen = 50, digits = 5)\r\n+\r\n+ if (!is.numeric(ids[[ivd]])) stop("Dependant variable is not numeric")\r\n+ if (!is.factor(ids[[ifixfact]])) stop("fixed factor is not a factor")\r\n+ if (!is.factor(ids[[itime]])) stop("Repeated factor is not a factor")\r\n+ if (!is.factor(ids[[isubject]])) stop("Random factor is not a factor")\r\n+\r\n+ ## factors\r\n+ time <- ids[[itime]]\r\n+ fixfact <- ids[[ifixfact]]\r\n+ subject <- ids[[isubject]]\r\n+ # dependant variables\r\n+ vd <- ids[[ivd]]\r\n+\r\n+ # argument of the function instead of re re-running ndim <- defColRes(ids, ifixfact, itime)\r\n+ # nfp : number of main factors + model infos (REML, varSubject) + normality test\r\n+ nfp <- ndim[1];\r\n+ # ncff number of comparison of the fixed factor\r\n+ nlff <- ndim[2]; ncff <- ndim[3]\r\n+ # nct number of comparison of the time factor\r\n+ nlt <- ndim[4]; nct <- ndim[5]\r\n+ # nci number of comparison of the interaction\r\n+ nli <- ndim[6]; nci <- ndim[7]\r\n+ # number of all lmer results\r\n+ nresT <- ncff + nct + nci\r\n+ ## initialization of the result vector (1 line)\r\n+ ## 4 * because nresf for : pvalues + Etimates + lower CI + Upper CI\r\n+ res <- data.frame(array(rep(NA, (nfp + 4 * nresT))))\r\n+ colnames(res)[1] <- "resultLM"\r\n+\r\n+ ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip\r\n+ ### after excluding NA, table function is used to seek subjects with only 1 data\r\n+ ids <- ids[!is.na(ids[[ivd]]), ]\r\n+ skip <- length(which(table(ids[[isubject]]) == 1))\r\n+\r\n+ if (skip == 0) {\r\n+\r\n+ mfl <- lmer(vd ~ time + fixfact + time:fixfact + (1 | subject), ids) # lmer remix\r\n+\r\n+ rsum <- summary(mfl, ddf = dffOption)\r\n+ ## test Shapiro Wilks on the residus of the model\r\n+ rShapiro <- shapiro.test(rsum$residuals)\r\n+ raov <- anova(mfl, ddf = dffOption)\r\n+ dlsm1 <- data.frame(difflsmeans(mfl, test.effs = NULL))\r\n+ ddlsm1 <- dlsm1\r\n+ ## save rownames and factor names\r\n+ rn <- rownames(ddlsm1)\r\n+ fn <- ddlsm1[, c(1, 2)]\r\n+ ## writing the results on a single line\r\n+ namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesFactPval <- paste("pvalue ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesInter <- rownames(ddlsm1)[-c(1:(nct + ncff))]\r\n+ namesEstimate <- paste("estimate ", namesInter)\r\n+ namespvalues <- paste("pvalue ", namesInter)\r\n+ namesFactprinc <- c("pval_time", "pval_trt", "pval_inter")\r\n+ namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+\r\n+ namesFactLowerCI <- paste("lowerCI ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesLowerCI <- paste("lowerCI ", namesInter, sep = "")\r\n+\r\n+ namesFactUpperCI <- paste("UpperCI ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesUpperCI <- paste("UpperCI ", namesInter, sep = "")\r\n+\r\n+\r\n+ ### lmer results on 1 vector row\r\n+ # pvalue of shapiro Wilks test of the residuals\r\n+ res[1, ] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals"\r\n+ res[2, ] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"\r\n+ res[3, ] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"\r\n+ '..b'utof), c((ndim[1] + ndim[5] + 1):(ndim[1] + ndim[5] + ndim[3]))] <- NA\r\n+ ## interaction effect\r\n+ resLM[which(resLM[, firstpval + 2] > pvalCutof), c((lastpval + ndim[5] + ndim[3] + 1):(lastpval + ndim[5] + ndim[3] + ndim[7]))] <- 0\r\n+ resLM[which(resLM[, firstpval + 2] > pvalCutof), c((ndim[1] + ndim[5] + ndim[3] + 1):(ndim[1] + ndim[5] + ndim[3] + ndim[7]))] <- NA\r\n+ } else {\r\n+ ## time effect only\r\n+ resLM[which(resLM[, firstpval] > pvalCutof), c((lastpval + 1):(lastpval + ndim[5]))] <- 0\r\n+ resLM[which(resLM[, firstpval] > pvalCutof), c((firstpval + 1):(firstpval + ndim[5]))] <- NA\r\n+ }\r\n+\r\n+ ## for each variable, estimates plots are performed if at least one factor is significant after p-value correction\r\n+ pdf(pdfE, onefile = TRUE, height = 15, width = 30)\r\n+\r\n+ ## for each variable (in row)\r\n+ for (i in seq_len(nrow(resLM))) {\r\n+\r\n+ ## if any fixed factor + time factor\r\n+ if (ifixfact > 0)\r\n+\r\n+ ## if any main factor after p-value correction is significant -> plot estimates and time course\r\n+ if (length(which(resLM[i, c(4:6)] < pvalCutof)) > 0) {\r\n+\r\n+ ## Plot of time course by fixfact : data prep with factors and quantitative var to be plot\r\n+ subv <- dslm[, colnames(dslm) == rownames(resLM)[i]]\r\n+ subds <- data.frame(dslm[[ifixfact]], dslm[[itime]], dslm[[isubject]], subv)\r\n+ libvar <- c(fixfact, time, subject)\r\n+ colnames(subds) <- c(libvar, rownames(resLM)[i])\r\n+\r\n+ ## Plot of estimates with error bars for all fixed factors and interaction\r\n+ rddlsm1 <- t(resLM[i, ])\r\n+ pval <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "pvalue"]\r\n+ esti <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "estima"]\r\n+ loci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "lowerC"]\r\n+ upci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "UpperC"]\r\n+ rddlsm1 <- data.frame(pval, esti, loci, upci, factorRow)\r\n+ colnames(rddlsm1) <- c("p.value", "Estimate", "Lower.CI", "Upper.CI", colnames(factorRow))\r\n+ rownames(rddlsm1) <- labelRow\r\n+\r\n+ ## function for plotting these 2 graphs\r\n+ plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)\r\n+\r\n+ }\r\n+\r\n+ ## if only a time factor\r\n+ if (ifixfact == 0)\r\n+\r\n+ ## if time factor after p-value correction is significant -> plot time course\r\n+ if (length(which(resLM[i, 4] < pvalCutof)) > 0) {\r\n+\r\n+ ## Plot of time course : data prep with factors and quantitative var to be plot\r\n+ subv <- dslm[, colnames(dslm) == rownames(resLM)[i]]\r\n+ subds <- data.frame(dslm[[itime]], dslm[[isubject]], subv)\r\n+ libvar <- c(time, subject)\r\n+ colnames(subds) <- c(libvar, rownames(resLM)[i])\r\n+\r\n+ ## Plot of estimates with error bars for all fixed factors and interaction\r\n+ rddlsm1 <- t(resLM[i, ])\r\n+ pval <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "pvalue"]\r\n+ esti <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "estima"]\r\n+ loci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "lowerC"]\r\n+ upci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "upperC"]\r\n+ rddlsm1 <- data.frame(pval, esti, loci, upci, factorRow)\r\n+ colnames(rddlsm1) <- c("p.value", "Estimate", "Lower.CI", "Upper.CI", colnames(factorRow))\r\n+ rownames(rddlsm1) <- labelRow\r\n+\r\n+ ## function for plotting these 2 graphs\r\n+ plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)\r\n+\r\n+ }\r\n+\r\n+ }\r\n+ dev.off()\r\n+\r\n+ ## return result file with pvalues and estimates (exclude confidence interval used for plotting)\r\n+ iCI <- which(substr(colnames(resLM), 4, 7) == "erCI")\r\n+ resLM <- resLM[, - iCI]\r\n+ resLM <- cbind(varids, resLM)\r\n+ return(resLM)\r\n+}\r\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/mixmodel_wrapper.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/mixmodel_wrapper.R Mon May 16 12:31:58 2022 +0000 |
| [ |
| @@ -0,0 +1,200 @@ +#!/usr/bin/env Rscript + +library(lme4) ## mixed model computing +library(Matrix) +library(MASS) +library(lmerTest) ## computing pvalue and lsmeans from results of lme4 package +library(multtest) ## multiple testing +library(ggplot2) +library(gridExtra) +library(grid) + +source_local <- function(fname) { + argv <- commandArgs(trailingOnly = FALSE) + base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) + source(paste(base_dir, fname, sep = "/")) +} + +source_local("mixmodel_script.R") +source_local("diagmfl.R") + +parse_args <- function() { + args <- commandArgs() + start <- which(args == "--args")[1] + 1 + if (is.na(start)) { + return(list()) + } + seq_by2 <- seq(start, length(args), by = 2) + result <- as.list(args[seq_by2 + 1]) + names(result) <- args[seq_by2] + return(result) +} + +argVc <- unlist(parse_args()) + +##------------------------------ +## Initializing +##------------------------------ + +## options +##-------- + +strAsFacL <- options()$stringsAsFactors +options(stringsAsFactors = FALSE) + +## constants +##---------- + +modNamC <- "mixmodel" ## module name + +topEnvC <- environment() +flagC <- "\n" + +## functions +##---------- + +flgF <- function(tesC, + envC = topEnvC, + txtC = NA) { ## management of warning and error messages + + tesL <- eval(parse(text = tesC), envir = envC) + + if (!tesL) { + + sink(NULL) + stpTxtC <- ifelse(is.na(txtC), + paste0(tesC, " is FALSE"), + txtC) + + stop(stpTxtC, + call. = FALSE) + + } + +} ## flgF + +## log file +##--------- + +sink(argVc["information"]) + +cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "") +cat("\nParameters used:\n\n") +print(argVc) +cat("\n\n") + +## loading +##-------- + +datMN <- t(as.matrix(read.table(argVc["dataMatrix_in"], + check.names = FALSE, + header = TRUE, + row.names = 1, + sep = "\t"))) + +samDF <- read.table(argVc["sampleMetadata_in"], + check.names = FALSE, + header = TRUE, + row.names = 1, + sep = "\t") + +varDF <- read.table(argVc["variableMetadata_in"], + check.names = FALSE, + header = TRUE, + row.names = 1, + sep = "\t") + + +## checking +##--------- + +flgF("identical(rownames(datMN), rownames(samDF))", txtC = "Column names of the dataMatrix are not identical to the row names of the sampleMetadata; check your data with the 'Check Format' module in the 'Quality Control' section") +flgF("identical(colnames(datMN), rownames(varDF))", txtC = "Row names of the dataMatrix are not identical to the row names of the variableMetadata; check your data with the 'Check Format' module in the 'Quality Control' section") + +flgF("argVc['time'] %in% colnames(samDF)", txtC = paste0("Required time factor '", argVc["time"], "' could not be found in the column names of the sampleMetadata")) +flgF("argVc['subject'] %in% colnames(samDF)", txtC = paste0("Required subject factor '", argVc["subject"], "' could not be found in the column names of the sampleMetadata")) + +flgF("mode(samDF[, argVc['time']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc["time"], "' column of the sampleMetadata should contain either number only, or character only")) +flgF("mode(samDF[, argVc['subject']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc["subject"], "' column of the sampleMetadata should contain either number only, or character only")) + +flgF("argVc['adjC'] %in% c('holm', 'hochberg', 'hommel', 'bonferroni', 'BH', 'BY', 'fdr', 'none')") +flgF("argVc['trf'] %in% c('none', 'log10', 'log2')") + +flgF("0 <= as.numeric(argVc['thrN']) && as.numeric(argVc['thrN']) <= 1", txtC = "(corrected) p-value threshold must be between 0 and 1") +flgF("argVc['diaR'] %in% c('no', 'yes')") + + +##------------------------------ +## Formating +##------------------------------ + +if (argVc["dff"] == "Satt") { + dffmeth <- "Satterthwaite" +} else { + dffmeth <- "Kenward-Roger" +} + + +##------------------------------ +## Computation +##------------------------------ + + +varDFout <- lmixedm(datMN = datMN, + samDF = samDF, + varDF = varDF, + fixfact = argVc["fixfact"], + time = argVc["time"], + subject = argVc["subject"], + logtr = argVc["trf"], + pvalCutof = argVc["thrN"], + pvalcorMeth = argVc["adjC"], + dffOption = dffmeth, + visu = argVc["diaR"], + least.confounded = FALSE, + outlier.limit = 3, + pdfC = argVc["out_graph_pdf"], + pdfE = argVc["out_estim_pdf"] + ) + + + + +##------------------------------ +## Rounding +##------------------------------ + +if (argVc["rounding"] == "yes") { + varDFout[, which(!(colnames(varDFout) %in% colnames(varDF)))] <- apply(varDFout[, which(!(colnames(varDFout) %in% colnames(varDF)))], 2, round, digits = as.numeric(argVc["decplaces"])) +} + +##------------------------------ +## Ending +##------------------------------ + + +## saving +##-------- + +varDFout <- cbind.data.frame(variableMetadata = rownames(varDFout), + varDFout) + +write.table(varDFout, + file = argVc["variableMetadata_out"], + quote = FALSE, + row.names = FALSE, + sep = "\t") + +## closing +##-------- + +cat("\n\nEnd of '", modNamC, "' Galaxy module call: ", + as.character(Sys.time()), "\n", sep = "") +cat("\nInformation about R (version, Operating System, attached or loaded packages):\n\n") +sessionInfo() + +sink() + +options(stringsAsFactors = strAsFacL) + +rm(list = ls()) |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/OneFactor_dataMatrix.tsv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/OneFactor_dataMatrix.tsv Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,23 @@ +variableMetadata s1D2 s2D2 s3D2 s1D3 s2D3 s3D3 s1D4 s2D4 s3D4 s1D5 s2D5 s3D5 s1D6 s2D6 s3D6 +V1 1.06E+06 1.07E+06 8.93E+05 9.11E+05 8.07E+05 6.90E+05 7.07E+05 7.66E+05 6.82E+05 6.25E+05 7.25E+05 6.78E+05 7.59E+05 6.70E+05 7.41E+05 +V2 4.67E+06 4.78E+06 4.17E+06 4.84E+06 4.50E+06 4.41E+06 6.96E+06 7.05E+06 6.65E+06 4.12E+06 4.42E+06 4.39E+06 1.26E+07 1.28E+07 1.32E+07 +V3 9.03E+05 8.44E+05 7.31E+05 9.97E+05 7.97E+05 7.47E+05 7.18E+05 8.28E+05 7.81E+05 9.08E+05 8.35E+05 8.37E+05 1.12E+06 1.11E+06 1.16E+06 +V4 1.62E+05 1.27E+05 1.21E+05 1.41E+05 8.48E+04 1.06E+05 9.72E+04 1.52E+05 1.27E+05 1.29E+05 1.42E+05 1.17E+05 1.72E+05 1.67E+05 1.87E+05 +V5 5.22E+04 6.20E+04 5.55E+04 4.78E+04 3.41E+04 5.84E+04 6.30E+04 1.11E+05 1.16E+05 4.59E+04 3.34E+04 4.08E+04 9.14E+04 6.78E+04 1.08E+05 +V6 1.68E+05 1.83E+05 2.76E+05 2.30E+05 1.62E+05 3.13E+05 4.47E+05 5.34E+05 6.59E+05 2.38E+05 2.30E+05 1.81E+05 4.93E+05 3.83E+05 6.30E+05 +V7 8.85E+06 8.39E+06 8.56E+06 5.51E+06 5.07E+06 4.90E+06 2.90E+06 3.10E+06 3.01E+06 2.36E+06 2.40E+06 2.30E+06 5.06E+06 4.94E+06 4.93E+06 +V8 2.53E+05 2.26E+05 1.90E+05 2.75E+05 2.67E+05 2.22E+05 2.12E+05 2.43E+05 1.72E+05 1.51E+05 1.56E+05 1.50E+05 2.49E+05 2.22E+05 2.48E+05 +V9 7.91E+05 7.16E+05 6.57E+05 6.08E+05 4.62E+05 4.88E+05 5.49E+05 6.13E+05 5.54E+05 7.56E+05 7.81E+05 7.59E+05 1.47E+06 1.46E+06 1.60E+06 +V10 1.76E+04 1.79E+04 1.86E+04 4.00E+04 4.27E+04 3.30E+04 5.34E+04 6.26E+04 6.09E+04 1.04E+05 8.24E+04 7.16E+04 2.09E+05 2.07E+05 2.19E+05 +V11 1.03E+05 9.28E+04 9.34E+04 8.70E+04 5.73E+04 8.05E+04 1.82E+06 1.98E+06 1.76E+06 8.41E+04 7.80E+04 7.31E+04 5.62E+04 5.01E+04 7.48E+04 +V12 4.48E+04 6.71E+04 4.27E+04 3.54E+04 3.04E+04 2.73E+04 2.47E+04 2.85E+04 3.26E+04 2.11E+04 1.55E+04 2.71E+04 3.44E+04 4.14E+04 5.87E+04 +V13 4.48E+04 6.71E+04 4.27E+04 3.54E+04 3.04E+04 2.73E+04 2.47E+04 2.85E+04 3.26E+04 2.11E+04 1.55E+04 2.71E+04 3.44E+04 4.14E+04 5.87E+04 +V14 1.16E+04 6.64E+03 1.29E+04 2.59E+04 1.55E+04 2.21E+04 3.15E+04 4.05E+04 3.10E+04 4.37E+04 2.90E+04 2.67E+04 6.62E+04 1.00E+05 8.75E+04 +V15 3.04E+05 3.05E+05 2.69E+05 2.57E+05 1.99E+05 2.16E+05 1.97E+05 1.97E+05 2.03E+05 2.65E+05 2.71E+05 2.60E+05 1.67E+05 1.89E+05 1.99E+05 +V16 2.69E+05 2.45E+05 2.18E+05 1.23E+05 1.27E+05 1.04E+05 2.91E+05 2.82E+05 2.79E+05 2.37E+05 2.11E+05 2.18E+05 5.60E+05 5.89E+05 6.08E+05 +V17 9.30E+04 8.55E+04 7.01E+04 4.84E+04 4.06E+04 3.14E+04 8.77E+04 1.03E+05 9.48E+04 9.88E+04 6.91E+04 7.81E+04 3.19E+05 3.25E+05 3.76E+05 +V18 2.30E+04 3.87E+04 2.01E+04 9.56E+03 1.27E+04 1.69E+04 2.38E+04 2.66E+04 2.99E+04 5.98E+04 4.93E+04 5.13E+04 1.08E+05 1.30E+05 1.36E+05 +V19 6.25E+04 7.43E+04 4.71E+04 1.62E+04 1.05E+04 1.57E+04 3.73E+04 3.08E+04 3.67E+04 3.50E+04 3.90E+04 2.78E+04 4.64E+05 1.60E+05 1.75E+05 +V20 1.72E+04 1.58E+04 1.35E+04 1.20E+04 1.02E+04 8.85E+03 1.35E+04 1.71E+04 1.62E+04 1.64E+04 8.41E+03 6.02E+03 5.07E+04 7.33E+04 6.57E+04 +V21 1.86E+06 1.90E+06 1.90E+06 2.05E+06 2.06E+06 2.07E+06 1.97E+06 1.94E+06 1.90E+06 1.64E+06 1.67E+06 1.78E+06 1.65E+06 1.73E+06 1.77E+06 +V22 1.52E+04 8.17E+03 1.37E+04 1.39E+04 1.34E+04 1.36E+04 1.40E+04 1.67E+04 1.60E+04 1.33E+04 1.69E+04 1.20E+04 1.05E+04 1.15E+04 1.11E+04 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/OneFactor_sampleMetadata.tsv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/OneFactor_sampleMetadata.tsv Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,16 @@ +sampleMetadata time subject +s1D2 D2 s1 +s2D2 D2 s2 +s3D2 D2 s3 +s1D3 D3 s1 +s2D3 D3 s2 +s3D3 D3 s3 +s1D4 D4 s1 +s2D4 D4 s2 +s3D4 D4 s3 +s1D5 D5 s1 +s2D5 D5 s2 +s3D5 D5 s3 +s1D6 D6 s1 +s2D6 D6 s2 +s3D6 D6 s3 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/OneFactor_variableMetadata.tsv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/OneFactor_variableMetadata.tsv Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,23 @@ +variableMetadata name +V1 name1 +V2 name2 +V3 name3 +V4 name4 +V5 name5 +V6 name6 +V7 name7 +V8 name8 +V9 name9 +V10 name10 +V11 name11 +V12 name12 +V13 name13 +V14 name14 +V15 name15 +V16 name16 +V17 name17 +V18 name18 +V19 name19 +V20 name20 +V21 name21 +V22 name22 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/TwoFactor_dataMatrix.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/TwoFactor_dataMatrix.txt Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,6 @@ +variable S01T0 S01T1 S01T2 S01T3 S02T0 S02T1 S02T2 S02T3 S03T0 S03T1 S03T2 S03T3 S04T0 S04T1 S04T2 S04T3 S05T0 S05T1 S05T2 S05T3 S06T0 S06T1 S06T2 S06T3 S07T0 S07T1 S07T2 S07T3 S08T0 S08T1 S08T2 S08T3 S09T0 S09T1 S09T2 S09T3 S10T0 S10T1 S10T2 S10T3 S11T0 S11T1 S11T2 S11T3 S12T0 S12T1 S12T2 S12T3 S13T0 S13T1 S13T2 S13T3 S14T0 S14T1 S14T2 S14T3 S15T0 S15T1 S15T2 S15T3 S16T0 S16T1 S16T2 S16T3 S17T0 S17T1 S17T2 S17T3 S18T0 S18T1 S18T2 S18T3 +v0.4425_385 126.39 234.11 162.68 115.29 172.44 171.53 180.31 114.94 224.14 195.09 189.45 122.79 212.02 243.48 188.42 138.91 241.59 192.95 214.74 98.69 270.03 176.79 208.73 78.04 247.82 176.51 175.78 98.94 260.63 181.93 178.06 111.58 218.95 176.24 200.29 112.77 140.23 157.38 124.62 114.34 147.23 165.29 125.85 128.41 126.35 139.75 127.95 74.86 124.92 155.12 109.4 86.88 145.57 142 135.29 83.08 116.95 142.28 113.43 118.26 80.81 146.03 130.52 103.4 79.58 168 109.35 88.19 138.12 166.61 115.05 83.05 +v0.4498_625 108.9 133.53 128.66 205.38 121.48 112.23 142.97 199.35 127.44 102.45 173.31 252.14 123.96 96.93 118.7 234.54 146.33 109.1 178.69 233.61 154.27 138.81 157.97 213.18 176.06 108.79 142.11 116.14 138.78 131.42 172.65 214.14 159.19 102.17 137.67 211.04 120.1 69.85 65.75 92.68 121.22 70.04 76.06 91.89 81.18 73.14 80.93 92.69 110.68 66.84 79.72 92.89 93.5 62.98 87.39 98.32 84.59 65.18 99.23 36.94 85.1 80.49 95.17 115.22 113.47 65.32 53.29 113.25 106.68 74.04 66.5 88.96 +v0.4711_463 124.1 166.84 156.18 82.57 158.4 190.09 152.89 100.76 155.75 174.85 153.14 112.76 171.54 150.37 163.49 116.06 198.78 161.53 247.61 122.94 167.38 186.12 194.59 93.05 208.24 150.15 140.97 55.62 200.32 138.59 145.99 89.44 192.49 131.44 145.89 93.38 102.94 67.11 83.99 86.86 110.09 98.61 77.91 100.75 136.32 92.56 86.79 75.62 144.39 65.27 110.24 94.26 136.48 77.62 95.43 60.3 135.98 87.54 65.37 92.04 112.09 58.9 92.35 72.92 132.84 97.12 80.71 79.01 144.75 85.78 73.61 83.68 +v0.4715_659 688.13 739.41 529.62 354.46 828.48 660.76 545.32 318.51 798.21 714.6 617.95 400.94 815.2 696.85 611.31 354.28 909.08 792.05 650.85 386.48 897.32 729.1 638.99 420.19 967.92 765.14 579.41 369.82 929.8 759.62 617.5 339.81 989.48 781.37 606.77 386.24 289.67 359.43 240.13 294.74 314.71 338.43 203.29 264.08 349.01 365.76 245.06 283.41 327.8 354.82 228.03 263.83 300.57 362.51 242.41 282.21 280.17 344.99 253.88 273.4 304.37 325.89 292.93 292.98 317.8 318.16 205.93 245.92 310.23 330.51 230.27 270.3 +v0.4766_757 117.8 110.31 114.56 99.99 144.39 155.33 114.47 122 147.16 132.82 129.78 94.86 128.26 161.38 119.12 96.73 177.26 134.94 149.39 134.12 163.22 146.12 148.26 108.68 159.39 140.75 119.11 65.83 170.45 131.79 163.58 104.63 167.33 112.88 142.58 85.38 80.53 84.3 42.33 62.04 73.74 97.78 58 79.33 91.7 82.39 69.71 83.58 77.86 65.56 74.11 82.9 90.53 75.27 80.46 81.88 92.77 52.94 67.99 84.21 91.27 69.54 77.13 85.27 93.18 83.68 57.01 65.09 92.12 82.63 70.9 75.53 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/TwoFactor_sampleMetadata.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/TwoFactor_sampleMetadata.txt Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,73 @@ +sample sujet phenotype time +S01T0 S01 WT T0 +S01T1 S01 WT T1 +S01T2 S01 WT T2 +S01T3 S01 WT T3 +S02T0 S02 WT T0 +S02T1 S02 WT T1 +S02T2 S02 WT T2 +S02T3 S02 WT T3 +S03T0 S03 WT T0 +S03T1 S03 WT T1 +S03T2 S03 WT T2 +S03T3 S03 WT T3 +S04T0 S04 WT T0 +S04T1 S04 WT T1 +S04T2 S04 WT T2 +S04T3 S04 WT T3 +S05T0 S05 WT T0 +S05T1 S05 WT T1 +S05T2 S05 WT T2 +S05T3 S05 WT T3 +S06T0 S06 WT T0 +S06T1 S06 WT T1 +S06T2 S06 WT T2 +S06T3 S06 WT T3 +S07T0 S07 WT T0 +S07T1 S07 WT T1 +S07T2 S07 WT T2 +S07T3 S07 WT T3 +S08T0 S08 WT T0 +S08T1 S08 WT T1 +S08T2 S08 WT T2 +S08T3 S08 WT T3 +S09T0 S09 WT T0 +S09T1 S09 WT T1 +S09T2 S09 WT T2 +S09T3 S09 WT T3 +S10T0 S10 MT T0 +S10T1 S10 MT T1 +S10T2 S10 MT T2 +S10T3 S10 MT T3 +S11T0 S11 MT T0 +S11T1 S11 MT T1 +S11T2 S11 MT T2 +S11T3 S11 MT T3 +S12T0 S12 MT T0 +S12T1 S12 MT T1 +S12T2 S12 MT T2 +S12T3 S12 MT T3 +S13T0 S13 MT T0 +S13T1 S13 MT T1 +S13T2 S13 MT T2 +S13T3 S13 MT T3 +S14T0 S14 MT T0 +S14T1 S14 MT T1 +S14T2 S14 MT T2 +S14T3 S14 MT T3 +S15T0 S15 MT T0 +S15T1 S15 MT T1 +S15T2 S15 MT T2 +S15T3 S15 MT T3 +S16T0 S16 MT T0 +S16T1 S16 MT T1 +S16T2 S16 MT T2 +S16T3 S16 MT T3 +S17T0 S17 MT T0 +S17T1 S17 MT T1 +S17T2 S17 MT T2 +S17T3 S17 MT T3 +S18T0 S18 MT T0 +S18T1 S18 MT T1 +S18T2 S18 MT T2 +S18T3 S18 MT T3 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/TwoFactor_variableMetadata.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/TwoFactor_variableMetadata.txt Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,6 @@ +variable num +v0.4425_385 1 +v0.4498_625 2 +v0.4711_463 3 +v0.4715_659 4 +v0.4766_757 5 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel4repeated_measures/test-data/mixmodel_TwoFactor_variableMetadata.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel4repeated_measures/test-data/mixmodel_TwoFactor_variableMetadata.txt Mon May 16 12:31:58 2022 +0000 |
| b |
| @@ -0,0 +1,6 @@ +variableMetadata num Shapiro.pvalue.residuals Subject.Variance REML pval_time pval_trt pval_inter pvalue.timeT0...timeT1 pvalue.timeT0...timeT2 pvalue.timeT0...timeT3 pvalue.timeT1...timeT2 pvalue.timeT1...timeT3 pvalue.timeT2...timeT3 pvalue.fixfactMT...fixfactWT pvalue..timeT0.fixfactMT...timeT1.fixfactMT pvalue..timeT0.fixfactMT...timeT2.fixfactMT pvalue..timeT0.fixfactMT...timeT3.fixfactMT pvalue..timeT0.fixfactMT...timeT0.fixfactWT pvalue..timeT0.fixfactMT...timeT1.fixfactWT pvalue..timeT0.fixfactMT...timeT2.fixfactWT pvalue..timeT0.fixfactMT...timeT3.fixfactWT pvalue..timeT1.fixfactMT...timeT2.fixfactMT pvalue..timeT1.fixfactMT...timeT3.fixfactMT pvalue..timeT1.fixfactMT...timeT0.fixfactWT pvalue..timeT1.fixfactMT...timeT1.fixfactWT pvalue..timeT1.fixfactMT...timeT2.fixfactWT pvalue..timeT1.fixfactMT...timeT3.fixfactWT pvalue..timeT2.fixfactMT...timeT3.fixfactMT pvalue..timeT2.fixfactMT...timeT0.fixfactWT pvalue..timeT2.fixfactMT...timeT1.fixfactWT pvalue..timeT2.fixfactMT...timeT2.fixfactWT pvalue..timeT2.fixfactMT...timeT3.fixfactWT pvalue..timeT3.fixfactMT...timeT0.fixfactWT pvalue..timeT3.fixfactMT...timeT1.fixfactWT pvalue..timeT3.fixfactMT...timeT2.fixfactWT pvalue..timeT3.fixfactMT...timeT3.fixfactWT pvalue..timeT0.fixfactWT...timeT1.fixfactWT pvalue..timeT0.fixfactWT...timeT2.fixfactWT pvalue..timeT0.fixfactWT...timeT3.fixfactWT pvalue..timeT1.fixfactWT...timeT2.fixfactWT pvalue..timeT1.fixfactWT...timeT3.fixfactWT pvalue..timeT2.fixfactWT...timeT3.fixfactWT estimate.timeT0...timeT1 estimate.timeT0...timeT2 estimate.timeT0...timeT3 estimate.timeT1...timeT2 estimate.timeT1...timeT3 estimate.timeT2...timeT3 estimate.fixfactMT...fixfactWT estimate..timeT0.fixfactMT...timeT1.fixfactMT estimate..timeT0.fixfactMT...timeT2.fixfactMT estimate..timeT0.fixfactMT...timeT3.fixfactMT estimate..timeT0.fixfactMT...timeT0.fixfactWT estimate..timeT0.fixfactMT...timeT1.fixfactWT estimate..timeT0.fixfactMT...timeT2.fixfactWT estimate..timeT0.fixfactMT...timeT3.fixfactWT estimate..timeT1.fixfactMT...timeT2.fixfactMT estimate..timeT1.fixfactMT...timeT3.fixfactMT estimate..timeT1.fixfactMT...timeT0.fixfactWT estimate..timeT1.fixfactMT...timeT1.fixfactWT estimate..timeT1.fixfactMT...timeT2.fixfactWT estimate..timeT1.fixfactMT...timeT3.fixfactWT estimate..timeT2.fixfactMT...timeT3.fixfactMT estimate..timeT2.fixfactMT...timeT0.fixfactWT estimate..timeT2.fixfactMT...timeT1.fixfactWT estimate..timeT2.fixfactMT...timeT2.fixfactWT estimate..timeT2.fixfactMT...timeT3.fixfactWT estimate..timeT3.fixfactMT...timeT0.fixfactWT estimate..timeT3.fixfactMT...timeT1.fixfactWT estimate..timeT3.fixfactMT...timeT2.fixfactWT estimate..timeT3.fixfactMT...timeT3.fixfactWT estimate..timeT0.fixfactWT...timeT1.fixfactWT estimate..timeT0.fixfactWT...timeT2.fixfactWT estimate..timeT0.fixfactWT...timeT3.fixfactWT estimate..timeT1.fixfactWT...timeT2.fixfactWT estimate..timeT1.fixfactWT...timeT3.fixfactWT estimate..timeT2.fixfactWT...timeT3.fixfactWT +v0.4425_385 1 0.004001 0 -140.931331 0 0 0.000264 0.191909 0.277993 0 0.018695 0 0 0 0.001835 0.869646 0.006728 0 0 0 0.238471 0.002992 0 0.000044 0.003836 0.009486 0.000036 0.004238 0 0 0 0.180279 0 0 0 0.112058 0.170546 0.091717 0 0.745427 0 0 -0.030835 0.02558 0.192247 0.056415 0.223081 0.166666 -0.148875 -0.107496 -0.005448 0.092614 -0.252536 -0.206709 -0.195927 0.039343 0.102048 0.200111 -0.14504 -0.099213 -0.088431 0.146839 0.098062 -0.247088 -0.201261 -0.19048 0.044791 -0.34515 -0.299323 -0.288542 -0.053272 0.045827 0.056608 0.291879 0.010781 0.246052 0.23527 +v0.4498_625 2 0.000002 0 -120.084702 0 0 0.000801 0.000045 0.135374 0.047823 0.005607 0 0.000776 0 0.000128 0.004174 0.135488 0.000855 0.170794 0.000053 0 0.272071 0.012599 0 0.000001 0 0 0.149351 0 0.000049 0 0 0.000005 0.00515 0 0 0.038427 0.408218 0.000047 0.00448 0 0.000769 0.120489 0.041611 -0.055515 -0.078879 -0.176005 -0.097126 -0.249448 0.158724 0.115616 0.058829 -0.136156 -0.053901 -0.16855 -0.306015 -0.043108 -0.099895 -0.29488 -0.212625 -0.327274 -0.46474 -0.056787 -0.251772 -0.169517 -0.284166 -0.421632 -0.194985 -0.11273 -0.227379 -0.364845 0.082255 -0.032394 -0.16986 -0.114649 -0.252114 -0.137466 +v0.4711_463 3 0.474695 0.000093 -135.57425 0 0 0.000006 0.000013 0.000128 0 0.488485 0.000043 0.000004 0 0 0.000004 0.000001 0.000295 0.006152 0.002433 0.000335 0.51173 0.755875 0 0 0 0.046325 0.729076 0 0 0 0.173445 0 0 0 0.089918 0.319884 0.50045 0 0.745698 0 0 0.11757 0.100696 0.22644 -0.016874 0.108869 0.125743 -0.194248 0.200783 0.178185 0.190094 -0.13207 -0.097712 -0.108863 0.130715 -0.022598 -0.010689 -0.332853 -0.298495 -0.309645 -0.070068 0.011909 -0.310255 -0.275897 -0.287048 -0.04747 -0.322164 -0.287806 -0.298957 -0.059379 0.034358 0.023208 0.262785 -0.01115 0.228427 0.239578 +v0.4715_659 4 0.623926 0.000349 -238.084202 0 0 0 0.219999 0 0 0 0 0 0 0.001773 0 0.000326 0 0 0 0.000025 0 0 0 0 0 0.070602 0.000027 0 0 0 0 0 0 0 0 0.000006 0 0 0 0 0 0.012039 0.138015 0.211119 0.125976 0.19908 0.073104 -0.326171 -0.045354 0.116577 0.053052 -0.44462 -0.375188 -0.285168 -0.075434 0.161931 0.098406 -0.399266 -0.329834 -0.239814 -0.03008 -0.063525 -0.561197 -0.491765 -0.401744 -0.192011 -0.497672 -0.42824 -0.33822 -0.128486 0.069432 0.159453 0.369186 0.09002 0.299754 0.209734 +v0.4766_757 5 0.005552 0.000365 -149.179654 0.00001 0 0.00028 0.015826 0.000092 0.000002 0.083272 0.004665 0.23638 0 0.064552 0.000169 0.102895 0 0 0 0.062636 0.033569 0.819707 0 0 0 0.000447 0.019489 0 0 0 0 0 0 0 0.000898 0.106307 0.055879 0 0.755242 0.000055 0.000152 0.052746 0.090042 0.115324 0.037296 0.062579 0.025283 -0.224878 0.056412 0.121655 0.049578 -0.240112 -0.191032 -0.181683 -0.059041 0.065243 -0.006834 -0.296524 -0.247444 -0.238095 -0.115453 -0.072077 -0.361767 -0.312687 -0.303338 -0.180696 -0.28969 -0.24061 -0.231262 -0.10862 0.04908 0.058428 0.181071 0.009349 0.131991 0.122642 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel_script.R --- a/mixmodel_script.R Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| b'@@ -1,456 +0,0 @@\n-####### R functions to perform linear mixed model for repeated measures \n-####### on a multi var dataset using 3 files as used in W4M\n-##############################################################################################################\n-lmRepeated2FF <- function(ids, ifixfact, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]], \n- pvalCutof=0.05,dffOption, visu , tit = "", least.confounded = FALSE, outlier.limit =3) \n- {\n- ### function to perform linear mixed model with 1 Fixed factor + Time + random factor subject\n- ### based on lmerTest package providing functions giving the same results as SAS proc mixed\n- options(scipen = 50, digits = 5)\n- \n- if (!is.numeric(ids[[ivd]])) {stop("Dependant variable is not numeric")}\n- if (!is.factor(ids[[ifixfact]])) {stop("fixed factor is not a factor")}\n- if (!is.factor(ids[[itime]])) {stop("Repeated factor is not a factor")}\n- if (!is.factor(ids[[isubject]])) {stop("Random factor is not a factor")}\n- # a ce stade, il faudrait pr?voir des tests sur la validit? du plan d\'exp?rience\n- \n- time <- ids[[itime]]\n- fixfact <- ids[[ifixfact]]\n- subject <- ids[[isubject]]\n- vd <- ids[[ivd]]\n- \n- # argument of the function instead of re re-running ndim <- defColRes(ids,ifixfact,itime)\n- # nfp : number of main factors + model infos (REML, varSubject) + normality test \n- nfp <- ndim[1];\n- # ncff number of comparison of the fixed factor\n- nlff <- ndim[2]; ncff <- ndim[3]\n- # nct number of comparison of the time factor\n- nlt <- ndim[4] ; nct <- ndim[5]\n- # nci number of comparison of the interaction\n- nli <- ndim[6]; nci <- ndim[7]\n- # number of all lmer results\n- nresT <- ncff+nct+nci\n- ## initialization of the result vector (1 line)\n- ## 4 * because nresf for : pvalues + Etimates + lower CI + Upper CI\n- res <- data.frame(array(rep(NA,(nfp + 4 * nresT))))\n- colnames(res)[1] <- "resultLM"\n- \n- ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip\n- ### after excluding NA, table function is used to seek subjects with only 1 data\n- ids <- ids[!is.na(ids[[ivd]]),]\n- skip <- length(which(table(ids[[isubject]])==1))\n- \n- if (skip==0) {\n- \n- mfl <- lmer( vd ~ time + fixfact + time:fixfact + (1| subject), ids) # lmer remix\n- \n- # ## NL add \n- # ### DEPLACE APRES CALCUL PVALUES AJUSTEES ET NE FAIRE QUE SI AU MOINS 1 FACTEUR SIGNIFICATIF\n- # if(visu) diagmflF(mfl, title = tit, least.confounded = least.confounded, outlier.limit = outlier.limit)\n- # ## end of NL add\n- \n- rsum <- summary(mfl,ddf = dffOption)\n- ## test Shapiro Wilks on the residus of the model \n- rShapiro <- shapiro.test(rsum$residuals)\n- raov <- anova(mfl,ddf = dffOption)\n- dlsm1 <- data.frame(difflsmeans(mfl,test.effs=NULL))\n- ddlsm1 <- dlsm1\n- ## save rownames and factor names\n- rn <- rownames(ddlsm1)\n- fn <- ddlsm1[,c(1,2)]\n- ## writing the results on a single line\n- namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")\n- namesFactPval <- paste("pvalue ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")\n- namesInter <- rownames(ddlsm1)[-c(1:(nct+ncff))]\n- #ncI <- nchar(namesInter)\n- namesEstimate <- paste("estimate ",namesInter)\n- namespvalues <- paste("pvalue ",namesInter)\n- namesFactprinc <- c("pval_time","pval_trt","pval_inter")\n- namesFactEstim <- paste("estimate ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")\n- \n- namesFactLowerCI <- paste("lowerCI ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")\n- namesLowerCI <- paste("lowerCI ",namesInter,sep="")\n- \n- namesFactUpperCI <- paste("UpperCI ",rownames(ddlsm1)[c(1:(nct+ncff))],sep="")\n- namesUpperCI <- paste("UpperCI ",namesInter,sep="")\n- \n- \n- ### lmer results on 1 vector row\n- # pvalue of shapiro Wilks '..b'of NL modif\n- resLM[i-firstvar+1,] <- reslmer[[1]]\n- }, error=function(e){cat("ERROR : ",conditionMessage(e), "\\n");})\n- if (i==firstvar) {\n- colnames(resLM) <- colnames(reslmer[[1]])\n- labelRow <- reslmer[[2]]\n- factorRow <- reslmer[[3]]\n- }\n- }\n- ## for debug : ifixfact=1;itime=2;isubject=3; ivd=4;tit = colnames(dslm)[i]; ids <- subds\n-\n- \n- ## NL add\n- if(visu == "yes") dev.off()\n- ## end of NL add\n- \n- ## pvalue correction with p.adjust library multtest\n- ## Possible methods of pvalue correction\n- AdjustMeth <- c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY","fdr","none")\n- if (length(which(pvalcorMeth == AdjustMeth))==0) pvalcorMeth <- "none"\n-\n- if (pvalcorMeth !="none") {\n- for (k in firstpval:lastpval){\n- resLM[[k]]=p.adjust(resLM[[k]], method=pvalcorMeth, n=dim(resLM[k])[[1]])\n- \n- }\n- }\n-\n- ## for each variables, set pvalues to NA and estimates = 0 when pvalue of factor > pvalCutof value define by user\n- if (ifixfact>0) {\n- ## time effect\n- resLM[which(resLM[,firstpval]> pvalCutof),c((lastpval+1):(lastpval+ndim[5]))] <- 0\n- resLM[which(resLM[,firstpval]> pvalCutof),c((ndim[1]+1):(ndim[1]+ndim[5]))] <- NA\n- ## treatment effect\n- resLM[which(resLM[,firstpval+1]> pvalCutof),c((lastpval+ndim[5]+1):(lastpval+ndim[5]+ndim[3]))] <- 0\n- resLM[which(resLM[,firstpval+1]> pvalCutof),c((ndim[1]+ndim[5]+1):(ndim[1]+ndim[5]+ndim[3]))] <- NA\n- ## interaction effect\n- resLM[which(resLM[,firstpval+2]> pvalCutof),c((lastpval+ndim[5]+ndim[3]+1):(lastpval+ndim[5]+ndim[3]+ndim[7]))] <- 0\n- resLM[which(resLM[,firstpval+2]> pvalCutof),c((ndim[1]+ndim[5]+ndim[3]+1):(ndim[1]+ndim[5]+ndim[3]+ndim[7]))] <- NA\n- } else {\n- ## time effect only\n- resLM[which(resLM[,firstpval]> pvalCutof),c((lastpval+1):(lastpval+ndim[5]))] <- 0\n- resLM[which(resLM[,firstpval]> pvalCutof),c((firstpval+1):(firstpval+ndim[5]))] <- NA\n- }\n- \n- ## for each variable, estimates plots are performed if at least one factor is significant after p-value correction\n- pdf(pdfE, onefile=TRUE, height = 15, width = 30)\n- #par(mfrow=c(2,2))\n- \n- ## for each variable (in row) \n- for (i in 1:nrow(resLM)) {\n- #cat("\\n",rownames(resLM)[i])\n- ## if any main factor after p-value correction is significant -> plot estimates and time course\n- if (length(which(resLM[i,c(4:6)]<pvalCutof))>0) {\n- \n- ## Plot of time course by fixfact : data prep with factors and quantitative var to be plot\n- subv <- dslm[,colnames(dslm)==rownames(resLM)[i]]\n- subds <- data.frame(dslm[[ifixfact]],dslm[[itime]], dslm[[isubject]],subv)\n- #colnames(subds) <- c(colnames(dslm)[ifixfact],colnames(dslm)[itime],colnames(dslm)[isubject],rownames(resLM)[i] <- rownames(resLM)[i] )\n- libvar <- c(fixfact,time,subject)\n- colnames(subds) <- c(libvar,rownames(resLM)[i])\n- \n- ## Plot of estimates with error bars for all fixed factors and interaction\n- rddlsm1 <- t(resLM[i,])\n- pval <- rddlsm1[substr(rownames(rddlsm1),1,6)=="pvalue"]\n- esti <- rddlsm1[substr(rownames(rddlsm1),1,6)=="estima"]\n- loci <- rddlsm1[substr(rownames(rddlsm1),1,6)=="lowerC"]\n- upci <- rddlsm1[substr(rownames(rddlsm1),1,6)=="UpperC"]\n- rddlsm1 <- data.frame(pval,esti,loci,upci,factorRow)\n- colnames(rddlsm1) <- c("p.value","Estimate","Lower.CI","Upper.CI",colnames(factorRow))\n- rownames(rddlsm1) <- labelRow\n-\n- ## function for plotting these 2 graphs\n- plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)\n- \n- }\n- }\n- dev.off()\n- \n- ## return result file with pvalues and estimates (exclude confidence interval used for plotting)\n- iCI <- which(substr(colnames(resLM),4,7)=="erCI")\n- resLM <- resLM[,-iCI] \n- resLM <- cbind(varids,resLM)\n- return(resLM)\n-}\n-\n-\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 mixmodel_wrapper.R --- a/mixmodel_wrapper.R Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| [ |
| @@ -1,164 +0,0 @@ -#!/usr/bin/env Rscript - -library(batch) ## parseCommandArgs - -library(lme4) ## mixed model computing -library(Matrix) -library(lmerTest) ## computing pvalue and lsmeans from results of lme4 package -library(multtest) ## multiple testing - - -source_local <- function(fname){ - argv <- commandArgs(trailingOnly = FALSE) - base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) - source(paste(base_dir, fname, sep="/")) -} - -#source_local("univariate_script.R") -source_local("mixmodel_script.R") -source_local("diagmfl.R") - -argVc <- unlist(parseCommandArgs(evaluate=FALSE)) - -##------------------------------ -## Initializing -##------------------------------ - -## options -##-------- - -strAsFacL <- options()$stringsAsFactors -options(stringsAsFactors = FALSE) - -## constants -##---------- - -modNamC <- "mixmodel" ## module name - -topEnvC <- environment() -flagC <- "\n" - -## functions -##---------- - -flgF <- function(tesC, - envC = topEnvC, - txtC = NA) { ## management of warning and error messages - - tesL <- eval(parse(text = tesC), envir = envC) - - if(!tesL) { - - sink(NULL) - stpTxtC <- ifelse(is.na(txtC), - paste0(tesC, " is FALSE"), - txtC) - - stop(stpTxtC, - call. = FALSE) - - } - -} ## flgF - -## log file -##--------- - -sink(argVc["information"]) - -cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep="") - -## loading -##-------- - -datMN <- t(as.matrix(read.table(argVc["dataMatrix_in"], - check.names = FALSE, - header = TRUE, - row.names = 1, - sep = "\t"))) - -samDF <- read.table(argVc["sampleMetadata_in"], - check.names = FALSE, - header = TRUE, - row.names = 1, - sep = "\t") - -varDF <- read.table(argVc["variableMetadata_in"], - check.names = FALSE, - header = TRUE, - row.names = 1, - sep = "\t") - - -## checking -##--------- - -flgF("identical(rownames(datMN), rownames(samDF))", txtC = "Column names of the dataMatrix are not identical to the row names of the sampleMetadata; check your data with the 'Check Format' module in the 'Quality Control' section") -flgF("identical(colnames(datMN), rownames(varDF))", txtC = "Row names of the dataMatrix are not identical to the row names of the variableMetadata; check your data with the 'Check Format' module in the 'Quality Control' section") - -flgF("argVc['fixfact'] %in% colnames(samDF)", txtC = paste0("Required fixed factor '" , argVc['fixfact'], "' could not be found in the column names of the sampleMetadata")) -flgF("argVc['time'] %in% colnames(samDF)", txtC = paste0("Required time factor '" , argVc['time'] , "' could not be found in the column names of the sampleMetadata")) -flgF("argVc['subject'] %in% colnames(samDF)", txtC = paste0("Required subject factor '", argVc['subject'], "' could not be found in the column names of the sampleMetadata")) - -flgF("mode(samDF[, argVc['fixfact']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc['fixfact'], "' column of the sampleMetadata should contain either number only, or character only")) -flgF("mode(samDF[, argVc['time']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc['time'] , "' column of the sampleMetadata should contain either number only, or character only")) -flgF("mode(samDF[, argVc['subject']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc['subject'], "' column of the sampleMetadata should contain either number only, or character only")) - -flgF("argVc['adjC'] %in% c('holm', 'hochberg', 'hommel', 'bonferroni', 'BH', 'BY', 'fdr', 'none')") -flgF("argVc['trf'] %in% c('none', 'log10', 'log2')") - -flgF("0 <= as.numeric(argVc['thrN']) && as.numeric(argVc['thrN']) <= 1", txtC = "(corrected) p-value threshold must be between 0 and 1") -flgF("argVc['diaR'] %in% c('no', 'yes')") - - -##------------------------------ -## Computation -##------------------------------ - - -varDF <- lmixedm(datMN = datMN, - samDF = samDF, - varDF = varDF, - fixfact = argVc["fixfact"], - time = argVc["time"], - subject = argVc["subject"], - logtr = argVc['trf'], - pvalCutof = argVc['thrN'], - pvalcorMeth= argVc["adjC"], - dffOption = "Satterthwaite", - visu = argVc["diaR"], - least.confounded = FALSE, - outlier.limit = 3, - pdfC = argVc["out_graph_pdf"], - pdfE = argVc["out_estim_pdf"] - ) - - -##------------------------------ -## Ending -##------------------------------ - - -## saving -##-------- - -varDF <- cbind.data.frame(variableMetadata = rownames(varDF), - varDF) - -write.table(varDF, - file = argVc["variableMetadata_out"], - quote = FALSE, - row.names = FALSE, - sep = "\t") - -## closing -##-------- - -cat("\nEnd of '", modNamC, "' Galaxy module call: ", - as.character(Sys.time()), "\n", sep = "") - -sink() - -options(stringsAsFactors = strAsFacL) - -rm(list = ls()) |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 test-data/demo1_Samples.txt --- a/test-data/demo1_Samples.txt Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,25 +0,0 @@ -idsample idsujet treatment time -s1AT1 s1 A T1 -s2AT1 s2 A T1 -s3AT1 s3 A T1 -s4AT1 s4 A T1 -s5BT1 s5 B T1 -s6BT1 s6 B T1 -s7BT1 s7 B T1 -s8BT1 s8 B T1 -s1AT2 s1 A T2 -s2AT2 s2 A T2 -s3AT2 s3 A T2 -s4AT2 s4 A T2 -s5BT2 s5 B T2 -s6BT2 s6 B T2 -s7BT2 s7 B T2 -s8BT2 s8 B T2 -s1AT3 s1 A T3 -s2AT3 s2 A T3 -s3AT3 s3 A T3 -s4AT3 s4 A T3 -s5BT3 s5 B T3 -s6BT3 s6 B T3 -s7BT3 s7 B T3 -s8BT3 s8 B T3 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 test-data/demo1_matrix.txt --- a/test-data/demo1_matrix.txt Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| b'@@ -1,41 +0,0 @@\n-var\ts1AT1\ts2AT1\ts3AT1\ts4AT1\ts5BT1\ts6BT1\ts7BT1\ts8BT1\ts1AT2\ts2AT2\ts3AT2\ts4AT2\ts5BT2\ts6BT2\ts7BT2\ts8BT2\ts1AT3\ts2AT3\ts3AT3\ts4AT3\ts5BT3\ts6BT3\ts7BT3\ts8BT3\n-var1\t2.944804979\t49.75484579\t25.89657977\t31.21956172\t645.7733673\t207.8199481\t993.1398444\t771.2810657\t33.8574053\t10.27450194\t42.27358423\t0.259314085\t727.698667\t975.6024581\t501.887448\t101.55269\t26.00199056\t42.81100563\t24.51317508\t4.306905516\t255.7178467\t844.9974405\t590.486975\t586.4456671\n-var2\t3.667579463\t305.5106362\t689.8621337\t960.0407436\t515.7640893\t925.1943751\t414.9315685\t285.7495409\t948.0050691\t821.9239079\t686.8984302\t254.7828577\t866.587021\t520.530133\t835.3730004\t966.3284721\t165.9697054\t726.3462446\t72.57517711\t370.4785168\t349.7503837\t368.2889788\t592.7939878\t502.4401631\n-var3\t172.8780668\t522.9435356\t60.99162384\t981.0661058\t838.0443428\t245.4660659\t320.52462\t402.9176733\tNA\tNA\t924.6445422\t109.0789884\tNA\t885.9453218\t821.7969874\t677.6235027\tNA\tNA\t573.2824613\t809.8096851\tNA\t989.3551375\t438.6697135\t99.02984484\n-var4\t476.6964656\t877.5340214\t873.7968183\t260.3147842\t48.1446747\t49.99827418\t15.0500212\t2.454515607\t655.4373272\t158.3886794\t533.7596174\t75.0229305\t21.36363224\t7.115041484\t6.141056371\t43.55349275\t287.4224888\t833.7583586\t60.74044966\t63.68551611\t47.30457162\t33.36978463\t42.5284983\t32.03159121\n-var5\t947.6929416\t453.2570661\t815.904467\t347.4032396\t24.53297914\t787.2949824\t241.7110583\t538.24695\t506.6958296\t266.0208725\t69.09085633\t5.974715233\t337.6606509\t259.2015248\t849.6838592\t87.20042666\t426.9964488\t698.8356796\t339.6737985\t156.684501\t86.32011891\t30.19974119\t978.5537816\t199.8844321\n-var6\t31.82849796\t44.69044713\t38.31707386\t18.51302603\t214.3392118\t712.6254117\t505.3932955\t707.9672445\t46.62706509\t43.77479739\t8.943157968\t39.26259778\t662.8065096\t808.9399247\t352.764633\t205.7874917\t3.766673404\t5.418740715\t37.59618013\t24.94170082\t516.9694331\t599.1699534\t951.1241327\t531.0029822\n-var7\t48.58915493\t9.737940091\t3.070545533\t13.24393267\t35.71840252\t38.19313066\t13.71876739\t37.84443501\t219.4206035\t529.5333088\t53.1798515\t720.0361187\t752.6900262\t884.0021377\t30.21174379\t639.3453975\t240.6978432\t287.5449343\t548.0445628\t269.850493\t742.9348096\t278.06316\t486.3017402\t805.5343043\n-var8\t680.6370973\t71.8828146\t981.710732\t221.6863381\t12.03643502\t18.67598209\t34.62782874\t47.00794697\t160.3053655\t566.2693544\t829.457417\t499.9228264\t2.015472805\t30.21909563\t11.46036299\t13.81338086\t227.484243\t99.45633839\t681.8675495\t973.652002\t41.87877455\t19.80605008\t43.40502059\t44.99090783\n-var9\t635.471832\t144.3056969\t334.1090387\t194.5727886\t405.5066688\t413.3236386\t237.899457\t996.7447354\t394.5598965\t519.4043987\t334.3332952\t866.3076014\t70.92464897\t689.4803042\t641.9404835\t453.5002735\t946.3136083\t557.380772\t725.8870537\t496.8042416\t251.7627825\t818.34199\t427.2466856\t193.6539535\n-var10\t255.5452475\t349.8598576\t617.453113\t707.2373216\t370.3060574\t469.7181814\t950.858553\t16.30868366\t806.6413776\t533.2787143\t6.329383006\t474.5446616\t462.0279273\t635.0181573\t651.8771692\t230.5368272\t222.9092267\t933.4693531\t548.8575506\t930.0440301\t265.5413981\t108.8758062\t68.08971196\t51.93824899\n-var11\t838.9985259\t745.5941668\t654.6952748\t235.1477422\t596.3175073\t844.6677532\t605.1399495\t432.9177807\t222.7906951\t873.972494\t822.0492133\t36.99986055\t537.8279339\t822.7241457\t5.419705705\t722.8959106\t881.2207647\t635.6996394\t677.2407238\t714.8925607\t389.7490318\t807.083624\t436.4182718\t458.9015398\n-var12\t892.7197353\t335.7591691\t536.4785811\t74.7070765\t576.9435273\t849.4743011\t401.1719511\t956.6287203\t40.56447695\t20.80947493\t1.140173277\t13.15588959\t30.28143383\t41.51846899\t25.38649045\t20.38415994\t793.7861762\t747.996181\t964.4964041\t786.1129855\t923.5151221\t394.6684835\t826.0571324\t213.7310158\n-var13\t499.4591238\t575.155097\t743.3397119\t27.23255924\t127.2970324\t640.1411168\t139.5822564\t642.8816939\t505.788513\t749.2321969\t691.3797399\t885.1521409\t60.46414235\t591.4272376\t94.30606951\t970.3055793\t755.4305079\t932.7703718\t653.9442423\t940.2452863\t661.9597309\t255.5623762\t72.89694711\t37.33840067\n-var14\t764.8635745\t186.441139\t403.7151201\t560.2602192\t833.8914349\t191.505535\t44.84334646\t467.5189404\t7'..b'327\t0.075892853\t3.35103701\t970.023551\t322.2403265\t1099.36841\t74.43502979\t1785.51484\t1459.638754\t990.5479732\t622.4609113\t109.2705883\t195.858416\t1059.939364\t1169.690758\t279.4441535\t511.1997851\n-var28\t433.1795515\t659.3852008\t304.342599\t2.420033492\t602.8680774\t1663.739933\t1022.569927\t342.8899009\t396.0198013\t788.9557203\t944.3908309\t314.6700035\t782.0471554\t309.1982655\t420.9024191\t592.3394811\t912.6321664\t414.073431\t482.1836798\t300.8122585\t689.3589282\t407.5987977\t1970.506461\t1274.174077\n-var29\t75.1274556\t667.9528257\t949.5909143\t332.3614352\t797.9255526\t1759.324563\t1697.283532\t1353.092656\t190.6264431\t16.51271938\t563.7405071\t674.3045958\t229.4363482\t1203.567204\t1431.917446\t854.3388249\t726.1664838\t236.2903037\t715.0134159\t201.2570507\t705.0542982\t1288.135354\t1282.10826\t742.1465211\n-var30\t372.3648963\t386.2979775\t484.8330806\t607.8727647\t844.1851676\t1564.698819\t225.4369875\t222.5468039\t675.5516492\t900.2904825\t321.0200632\t124.8706177\t219.1452961\t1620.109962\t1443.686531\t30.20502128\t711.3623969\t472.7569176\t131.2124078\t4.11041311\t314.6987183\t1480.330689\t1335.262401\t1975.43864\n-var31\t830.7603594\t877.7867945\t167.2341465\t890.4469672\t475.8462996\t1286.985999\t325.0642312\t1206.54618\t871.5761917\t924.3084225\t375.9059868\t841.7447645\t223.2287944\t489.4455531\t1276.794717\t555.4166789\t665.9058972\t900.4792176\t606.8139413\t662.8645656\t1453.245611\t386.6570783\t28.38822884\t850.4592565\n-var32\t287.6665566\t704.2399753\t202.4551918\t749.5138905\t251.7461375\t27.86278674\t1587.838351\t725.2323834\t927.9858767\t879.604207\t533.5445114\t684.1694546\t926.6281568\t164.3194326\t226.4151943\t1448.842914\t313.5647655\t54.62111634\t772.0807022\t659.3092555\t1303.021777\t1946.40263\t250.5414262\t1999.222444\n-var33\t484.1877297\t242.788894\t391.1647338\t296.6169619\t566.7181555\t556.7362445\t1705.383419\t118.3629431\t196.6473393\t943.5570103\t514.8669868\t994.6891166\t1630.677314\t1743.953602\t1686.234573\t1025.512791\t139.1387544\t102.1118108\t138.6877081\t679.5824526\t1399.967774\t1964.658448\t1942.234917\t345.9381031\n-var34\t654.4772218\t855.843767\t710.7433264\t332.3445897\t822.7186717\t1732.030543\t1071.843154\t977.4774333\t555.832087\t947.4166002\t845.8471497\t52.26234852\t1759.519813\t252.0201757\t132.8152855\t460.9241952\t589.2474936\t793.5101148\t323.7116484\t698.2504426\t949.3640103\t549.1089047\t1268.608356\t1981.901916\n-var35\t649.8157994\t596.3701266\t459.9072448\t636.8172769\t756.8655665\t1682.670261\t1356.109792\t1582.682511\t984.4088071\t737.4489045\t431.9834858\t175.3068228\t1640.146811\t611.9998332\t1107.834995\t983.1653253\t636.2798018\t0.545485104\t498.2985486\t657.9607458\t1574.724131\t725.3778129\t1287.929168\t396.8659763\n-var36\t524.2048663\t329.6938613\t884.8469578\t890.8421732\t27.18495267\t914.3060082\t196.5344281\t625.4772393\t818.8061525\t722.1835868\t872.4183559\t972.937521\t762.369294\t793.0758578\t1222.96799\t360.5858065\t95.91319149\t654.6684817\t373.220399\t804.0331353\t1402.29363\t571.7247152\t1489.69958\t774.8183349\n-var37\t678.7777667\t363.5061746\t834.5242973\t684.8352169\t760.735293\t222.6371125\t60.32352661\t867.7113126\t479.1262946\t325.0865937\t434.5733939\t893.4304109\t1606.009688\t509.1794575\t1282.08164\t1094.820467\t12.68285865\t2.822794117\t940.643551\t315.5426037\t682.6764091\t1041.812735\t177.9454557\t848.1309378\n-var38\t706.0233496\t934.3787643\t820.7876871\t72.28857763\t493.9847378\t247.5833014\t638.2539566\t1344.694015\t688.708953\t343.7132988\t491.9668005\t757.1955777\t573.0037148\t1605.206211\t298.9895983\t1269.280218\t994.2564165\t207.6765578\t299.4243272\t956.917879\t771.6799955\t223.0914032\t1996.397524\t67.29379012\n-var39\t490.8321846\t438.2277554\t656.9503593\t933.3194445\t868.3614729\t366.8485261\t1513.524538\t764.5321233\t871.1015601\t569.8457348\t323.2530982\t484.6563769\t368.7677212\t833.4574507\t1390.149132\t78.84558674\t483.1216159\t309.9850905\t657.4770171\t587.362726\t308.5230175\t1303.562509\t1352.843606\t474.7223959\n-var40\t50.90655113\t69.18390457\t457.9171462\t405.0173308\t146.918983\t1142.527384\t1412.337297\t850.7689924\t35.11405349\t23.21297831\t792.7634193\t387.1686795\t838.6976213\t1603.972482\t1325.575054\t1448.659011\t31.93235194\t135.4397238\t75.06461922\t543.702984\t216.6262105\t220.2400455\t1875.783019\t1121.44313\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 test-data/demo1_variables.txt --- a/test-data/demo1_variables.txt Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,41 +0,0 @@ -var withNA -var1 0 -var2 0 -var3 1 -var4 0 -var5 0 -var6 0 -var7 0 -var8 0 -var9 0 -var10 0 -var11 0 -var12 0 -var13 0 -var14 0 -var15 0 -var16 0 -var17 0 -var18 0 -var19 0 -var20 0 -var21 0 -var22 0 -var23 0 -var24 0 -var25 0 -var26 0 -var27 0 -var28 0 -var29 0 -var30 0 -var31 0 -var32 0 -var33 0 -var34 0 -var35 0 -var36 0 -var37 0 -var38 0 -var39 0 -var40 0 |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 test-data/mixmodel_demo1_variables.txt --- a/test-data/mixmodel_demo1_variables.txt Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| b'@@ -1,41 +0,0 @@\n-variableMetadata\twithNA\tShapiro.pvalue.residuals\tSubject.Variance\tREML\tpval_time\tpval_trt\tpval_inter\tpvalue.timeT1...timeT2\tpvalue.timeT1...timeT3\tpvalue.timeT2...timeT3\tpvalue.fixfactA...fixfactB\tpvalue..timeT1.fixfactA...timeT2.fixfactA\tpvalue..timeT1.fixfactA...timeT3.fixfactA\tpvalue..timeT1.fixfactA...timeT1.fixfactB\tpvalue..timeT1.fixfactA...timeT2.fixfactB\tpvalue..timeT1.fixfactA...timeT3.fixfactB\tpvalue..timeT2.fixfactA...timeT3.fixfactA\tpvalue..timeT2.fixfactA...timeT1.fixfactB\tpvalue..timeT2.fixfactA...timeT2.fixfactB\tpvalue..timeT2.fixfactA...timeT3.fixfactB\tpvalue..timeT3.fixfactA...timeT1.fixfactB\tpvalue..timeT3.fixfactA...timeT2.fixfactB\tpvalue..timeT3.fixfactA...timeT3.fixfactB\tpvalue..timeT1.fixfactB...timeT2.fixfactB\tpvalue..timeT1.fixfactB...timeT3.fixfactB\tpvalue..timeT2.fixfactB...timeT3.fixfactB\testimate.timeT1...timeT2\testimate.timeT1...timeT3\testimate.timeT2...timeT3\testimate.fixfactA...fixfactB\testimate..timeT1.fixfactA...timeT2.fixfactA\testimate..timeT1.fixfactA...timeT3.fixfactA\testimate..timeT1.fixfactA...timeT1.fixfactB\testimate..timeT1.fixfactA...timeT2.fixfactB\testimate..timeT1.fixfactA...timeT3.fixfactB\testimate..timeT2.fixfactA...timeT3.fixfactA\testimate..timeT2.fixfactA...timeT1.fixfactB\testimate..timeT2.fixfactA...timeT2.fixfactB\testimate..timeT2.fixfactA...timeT3.fixfactB\testimate..timeT3.fixfactA...timeT1.fixfactB\testimate..timeT3.fixfactA...timeT2.fixfactB\testimate..timeT3.fixfactA...timeT3.fixfactB\testimate..timeT1.fixfactB...timeT2.fixfactB\testimate..timeT1.fixfactB...timeT3.fixfactB\testimate..timeT2.fixfactB...timeT3.fixfactB\n-var1\t0\t0.00152990892918407\t0\t254.537609570589\t0.989895720028428\t0.000138446468356235\t0.999849571595415\tNA\tNA\tNA\t0.000138448022388127\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t-575.690811941667\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var2\t0\t0.960464619805585\t0\t262.571810056621\t0.448325032764408\t0.602586916593328\t1\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var3\t1\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var4\t0\t0.0441977246388017\t7213.02082711591\t254.223292186154\t0.983047523518005\t0.0690465434631971\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var5\t0\t0.064287156921339\t17220.7805833911\t266.097162357616\t0.983047523518005\t0.850320276396892\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var6\t0\t0.0381978835199867\t0\t244.361560756404\t0.983047523518005\t0.0000168443985191467\t0.877942440820325\tNA\tNA\tNA\t0.0000168442499891452\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t-535.434188801917\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var7\t0\t0.193737436286724\t0\t254.906518619133\t0.0156061251704555\t0.293289541312568\t0.877942440820325\t0.0118965177919926\t0.0269815160377476\t0.977421297378864\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t-453.53785986075\t-432.3569423245\t21.18091753625\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var8\t0\t0.014695144502403\t15197.0828494013\t259.222125604721\t1\t0.0697161146068648\t1\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var9\t0\t0.148540192056018\t0.00000000018132043075478\t259.983587153763\t0.983047523518005\t0.818083737709546\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var10\t0\t0.768501224171673\t4386.50412060244\t262.180681429873\t0.983047523518005\t0.372676427326723\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var11\t0\t0.711106348589254\t19233.8026294981\t260.769418500537\t0.983047523518005\t0.819711630041906\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var12\t0\t0.388310543095299\t0\t254.652163866062\t0.000464327919921784\t1\t0.86214265019290'..b'51387312325\t0.983047523518005\t0.71138248160657\t0.924943894585718\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var22\t0\t0.607405931022872\t0\t263.535518930099\t0.983047523518005\t0.569060445862623\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var23\t0\t0.212615804058253\t0\t263.751117472845\t0.983047523518005\t0.644748172383697\t0.924943894585718\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var24\t0\t0.453100802782389\t29510.6568674354\t260.71769299974\t0.989895720028428\t0.853217576536927\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var25\t0\t0.511477701192744\t0\t286.65672225126\t0.983047523518005\t0.286161368832769\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var26\t0\t0.815661522629652\t0\t281.509414777202\t0.983047523518005\t0.0690465434631971\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var27\t0\t0.653920548984792\t0\t280.838740914533\t0.989895720028428\t0.150203459443103\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var28\t0\t0.719759503524831\t0\t277.41708583405\t0.983047523518005\t0.1969140094154\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var29\t0\t0.725675776865143\t86615.1480315815\t268.196774164521\t0.448325032764408\t0.121472365285413\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var30\t0\t0.857478473660295\t37081.0431925278\t286.162827544403\t0.983047523518005\t0.253534413626269\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var31\t0\t0.751744017902084\t0\t276.362228496729\t0.989895720028428\t1\t0.924943894585718\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var32\t0\t0.978108618401383\t0\t285.72256677923\t0.983047523518005\t0.306152094558791\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var33\t0\t0.99993316514364\t124822.847216428\t276.589351488647\t0.144068402464455\t0.121472365285413\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var34\t0\t0.242528626574342\t0\t281.137098865892\t0.983047523518005\t0.1969140094154\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var35\t0\t0.508321244841884\t0\t273.186257975745\t0.983047523518005\t0.0141111367336112\t0.924943894585718\tNA\tNA\tNA\t0.0141114522811781\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t-603.435761166334\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var36\t0\t0.316164521545131\t0.0000000010285854101153\t268.90084474405\t0.983047523518005\t0.644748172383697\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var37\t0\t0.953719843399109\t20144.5877011241\t271.448160048801\t0.806157065268918\t0.356762591033804\t0.862142650192906\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var38\t0\t0.428134779270157\t0\t285.508280238472\t0.989895720028428\t0.592583606736207\t0.924943894585718\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var39\t0\t0.868929253676507\t74721.1043100503\t271.9525623048\t0.983047523518005\t0.569060445862623\t0.924943894585718\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n-var40\t0\t0.967902408294312\t97407.8193511745\t276.870664881354\t0.983047523518005\t0.121472365285413\t0.877942440820325\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\t0\n' |
| b |
| diff -r 1422de181204 -r a3147e3d66e2 test-data/test_raw_variable_out.txt --- a/test-data/test_raw_variable_out.txt Wed Oct 10 05:18:42 2018 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| b'@@ -1,41 +0,0 @@\n-variableMetadata\tvar\twithNA\tShapiro.pvalue.residuals\tSubject.Variance\tREML\tpval_time\tpval_trt\tpval_inter\tpvalue.time.T1...T2\tpvalue.time.T1...T3\tpvalue.time.T2...T3\tpvalue.fixfact.A...B\tpvalue.T1.A....T2.A\tpvalue.T1.A....T3.A\tpvalue.T1.A....T1.B\tpvalue.T1.A....T2.B\tpvalue.T1.A....T3.B\tpvalue.T2.A....T3.A\tpvalue.T2.A....T1.B\tpvalue.T2.A....T2.B\tpvalue.T2.A....T3.B\tpvalue.T3.A....T1.B\tpvalue.T3.A....T2.B\tpvalue.T3.A....T3.B\tpvalue.T1.B....T2.B\tpvalue.T1.B....T3.B\tpvalue.T2.B....T3.B\testimate.time.T1...T2\testimate.time.T1...T3\testimate.time.T2...T3\testimate.fixfact.A...B\testimate.T1.A....T2.A\testimate.T1.A....T3.A\testimate.T1.A....T1.B\testimate.T1.A....T2.B\testimate.T1.A....T3.B\testimate.T2.A....T3.A\testimate.T2.A....T1.B\testimate.T2.A....T2.B\testimate.T2.A....T3.B\testimate.T3.A....T1.B\testimate.T3.A....T2.B\testimate.T3.A....T3.B\testimate.T1.B....T2.B\testimate.T1.B....T3.B\testimate.T2.B....T3.B\n-var1\tvar1\t0\t0.00152990892918407\t0\t254.537609570589\t0.989895081688929\t0.000142332953130264\t0.999851637979215\tNA\tNA\tNA\t0\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t-575.6908\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var2\tvar2\t0\t0.960464619805585\t0\t262.571810056621\t0.443498818238816\t0.602652253312485\t1\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var3\tvar3\t1\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var4\tvar4\t0\t0.0441977246388017\t7213.02082711591\t254.223292186154\t0.983046491866558\t0.0690481436570156\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var5\tvar5\t0\t0.064287156921339\t17220.7805833911\t266.097162357616\t0.983046491866558\t0.850330775367304\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var6\tvar6\t0\t0.0381978835199867\t0\t244.361560756404\t0.983046491866558\t1.69724227738755e-05\t0.877939640795251\tNA\tNA\tNA\t0\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t-535.4342\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var7\tvar7\t0\t0.193737436286724\t0\t254.906518619133\t0.0157977802240783\t0.293529919182109\t0.877939640795251\t0.012\t0.028\t0.977485714285714\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t-453.5379\t-432.3569\t21.1809\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var8\tvar8\t0\t0.014695144502403\t15197.0828494013\t259.222125604721\t1\t0.0697160454979429\t1\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var9\tvar9\t0\t0.148540192056018\t1.8132043075478e-10\t259.983587153763\t0.983046491866558\t0.818040560247076\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var10\tvar10\t0\t0.768501224171673\t4386.50412060244\t262.180681429873\t0.983046491866558\t0.372670267631705\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var11\tvar11\t0\t0.711106348589254\t19233.8026294981\t260.769418500537\t0.983046491866558\t0.819708219928867\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var12\tvar12\t0\t0.388310543095299\t0\t254.652163866062\t0.000451504869487529\t1\t0.861541599913403\t0.002\t0.908\t0\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t553.8303\t-128.3101\t-682.1404\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var13\tvar13\t0\t0.846722054364422\t0\t263.195387302844\t0.983046491866558\t0.102587754537138\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var14\tvar14\t0\t0.0360856422369758\t0\t255.35073473538\t0.0804913814946417\t1\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tN'..b'44\t0.924958848665671\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var24\tvar24\t0\t0.453100802782389\t29510.6568674354\t260.71769299974\t0.989895081688929\t0.853232770669574\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var25\tvar25\t0\t0.511477701192744\t0\t286.65672225126\t0.983046491866558\t0.285791545615653\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var26\tvar26\t0\t0.815661522629652\t0\t281.509414777202\t0.983046491866558\t0.0688280391920837\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var27\tvar27\t0\t0.653920548984792\t0\t280.838740914533\t0.989895081688929\t0.14979041564372\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var28\tvar28\t0\t0.719759503524831\t0\t277.41708583405\t0.983046491866558\t0.196622531113559\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var29\tvar29\t0\t0.725675776865143\t86615.1480315815\t268.196774164521\t0.443498818238816\t0.121470044896006\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var30\tvar30\t0\t0.857478473660295\t37081.0431925278\t286.162827544403\t0.983046491866558\t0.253530786238282\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var31\tvar31\t0\t0.751744017902084\t0\t276.362228496729\t0.989895081688929\t1\t0.924958848665671\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var32\tvar32\t0\t0.978108618401383\t0\t285.72256677923\t0.983046491866558\t0.306024424517496\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var33\tvar33\t0\t0.99993316514364\t124822.847216428\t276.589351488647\t0.147948925847613\t0.121470044896006\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var34\tvar34\t0\t0.242528626574342\t0\t281.137098865892\t0.983046491866558\t0.196622531113559\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var35\tvar35\t0\t0.508321244841884\t0\t273.186257975745\t0.983046491866558\t0.0141133867640431\t0.924958848665671\tNA\tNA\tNA\t0.0146666666666667\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\t-603.4358\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var36\tvar36\t0\t0.316164521545126\t2.61433151397349e-10\t268.90084474405\t0.983046491866558\t0.644803266722344\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var37\tvar37\t0\t0.953719843399109\t20144.5877011241\t271.448160048801\t0.800264575774149\t0.356750278198418\t0.861541599913403\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var38\tvar38\t0\t0.428134779270157\t0\t285.508280238472\t0.989895081688929\t0.592722965710858\t0.924958848665671\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var39\tvar39\t0\t0.868929253676507\t74721.1043100503\t271.9525623048\t0.983046491866558\t0.568870942400682\t0.924958848665671\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n-var40\tvar40\t0\t0.967902408294312\t97407.8193511745\t276.870664881354\t0.983046491866558\t0.121470044896006\t0.877939640795251\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\tNA\n' |