Mercurial > repos > george-weingart > maaslin
annotate maaslin-4450aa4ecc84/src/lib/Misc.R @ 1:a87d5a5f2776
Uploaded the version running on the prod server
author | george-weingart |
---|---|
date | Sun, 08 Feb 2015 23:08:38 -0500 |
parents | |
children |
rev | line source |
---|---|
1
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
1 ##################################################################################### |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
2 #Copyright (C) <2012> |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
3 # |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
4 #Permission is hereby granted, free of charge, to any person obtaining a copy of |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
5 #this software and associated documentation files (the "Software"), to deal in the |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
6 #Software without restriction, including without limitation the rights to use, copy, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
7 #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
8 #and to permit persons to whom the Software is furnished to do so, subject to |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
9 #the following conditions: |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
10 # |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
11 #The above copyright notice and this permission notice shall be included in all copies |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
12 #or substantial portions of the Software. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
13 # |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
14 #THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
15 #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
16 #PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
17 #HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
18 #OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
19 #SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
20 # |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
21 # This file is a component of the MaAsLin (Multivariate Associations Using Linear Models), |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
22 # authored by the Huttenhower lab at the Harvard School of Public Health |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
23 # (contact Timothy Tickle, ttickle@hsph.harvard.edu). |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
24 ##################################################################################### |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
25 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
26 ### Modified Code |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
27 ### This code is from the package agricolae by Felipe de Mendiburu |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
28 ### Modifications here are minimal and allow one to use the p.values from the post hoc comparisons |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
29 ### Authors do not claim credit for this solution only needed to modify code to use the output. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
30 kruskal <- function (y, trt, alpha = 0.05, p.adj = c("none", "holm", "hochberg", |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
31 "bonferroni", "BH", "BY", "fdr"), group = TRUE, main = NULL) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
32 { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
33 dfComparisons=NULL |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
34 dfMeans=NULL |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
35 dntStudent=NULL |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
36 dLSD=NULL |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
37 dHMean=NULL |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
38 name.y <- paste(deparse(substitute(y))) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
39 name.t <- paste(deparse(substitute(trt))) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
40 p.adj <- match.arg(p.adj) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
41 junto <- subset(data.frame(y, trt), is.na(y) == FALSE) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
42 N <- nrow(junto) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
43 junto[, 1] <- rank(junto[, 1]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
44 means <- tapply.stat(junto[, 1], junto[, 2], stat = "sum") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
45 sds <- tapply.stat(junto[, 1], junto[, 2], stat = "sd") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
46 nn <- tapply.stat(junto[, 1], junto[, 2], stat = "length") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
47 means <- data.frame(means, replication = nn[, 2]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
48 names(means)[1:2] <- c(name.t, name.y) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
49 ntr <- nrow(means) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
50 nk <- choose(ntr, 2) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
51 DFerror <- N - ntr |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
52 rs <- 0 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
53 U <- 0 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
54 for (i in 1:ntr) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
55 rs <- rs + means[i, 2]^2/means[i, 3] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
56 U <- U + 1/means[i, 3] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
57 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
58 S <- (sum(junto[, 1]^2) - (N * (N + 1)^2)/4)/(N - 1) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
59 H <- (rs - (N * (N + 1)^2)/4)/S |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
60 # cat("\nStudy:", main) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
61 # cat("\nKruskal-Wallis test's\nTies or no Ties\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
62 # cat("\nValue:", H) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
63 # cat("\ndegrees of freedom:", ntr - 1) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
64 p.chisq <- 1 - pchisq(H, ntr - 1) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
65 # cat("\nPvalue chisq :", p.chisq, "\n\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
66 DFerror <- N - ntr |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
67 Tprob <- qt(1 - alpha/2, DFerror) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
68 MSerror <- S * ((N - 1 - H)/(N - ntr)) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
69 means[, 2] <- means[, 2]/means[, 3] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
70 # cat(paste(name.t, ",", sep = ""), " means of the ranks\n\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
71 dfMeans=data.frame(row.names = means[, 1], means[, -1]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
72 if (p.adj != "none") { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
73 # cat("\nP value adjustment method:", p.adj) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
74 a <- 1e-06 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
75 b <- 1 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
76 for (i in 1:100) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
77 x <- (b + a)/2 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
78 xr <- rep(x, nk) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
79 d <- p.adjust(xr, p.adj)[1] - alpha |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
80 ar <- rep(a, nk) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
81 fa <- p.adjust(ar, p.adj)[1] - alpha |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
82 if (d * fa < 0) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
83 b <- x |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
84 if (d * fa > 0) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
85 a <- x |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
86 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
87 Tprob <- qt(1 - x/2, DFerror) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
88 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
89 nr <- unique(means[, 3]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
90 if (group) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
91 Tprob <- qt(1 - alpha/2, DFerror) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
92 # cat("\nt-Student:", Tprob) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
93 # cat("\nAlpha :", alpha) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
94 dntStudent=Tprob |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
95 dAlpha=alpha |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
96 if (length(nr) == 1) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
97 LSD <- Tprob * sqrt(2 * MSerror/nr) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
98 # cat("\nLSD :", LSD, "\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
99 dLSD=LSD |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
100 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
101 else { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
102 nr1 <- 1/mean(1/nn[, 2]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
103 LSD1 <- Tprob * sqrt(2 * MSerror/nr1) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
104 # cat("\nLSD :", LSD1, "\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
105 dLSD =LSD1 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
106 # cat("\nHarmonic Mean of Cell Sizes ", nr1) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
107 dHMean=nr1 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
108 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
109 # cat("\nMeans with the same letter are not significantly different\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
110 # cat("\nGroups, Treatments and mean of the ranks\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
111 output <- order.group(means[, 1], means[, 2], means[, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
112 3], MSerror, Tprob, std.err = sqrt(MSerror/means[, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
113 3])) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
114 dfComparisons=order.group(means[, 1], means[, 2], means[, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
115 3], MSerror, Tprob, std.err = sqrt(MSerror/means[, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
116 3])) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
117 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
118 if (!group) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
119 comb <- combn(ntr, 2) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
120 nn <- ncol(comb) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
121 dif <- rep(0, nn) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
122 LCL <- dif |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
123 UCL <- dif |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
124 pvalue <- dif |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
125 sdtdif <- dif |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
126 for (k in 1:nn) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
127 i <- comb[1, k] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
128 j <- comb[2, k] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
129 if (means[i, 2] < means[j, 2]) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
130 comb[1, k] <- j |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
131 comb[2, k] <- i |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
132 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
133 dif[k] <- abs(means[i, 2] - means[j, 2]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
134 sdtdif[k] <- sqrt(S * ((N - 1 - H)/(N - ntr)) * (1/means[i, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
135 3] + 1/means[j, 3])) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
136 pvalue[k] <- 2 * round(1 - pt(dif[k]/sdtdif[k], DFerror), |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
137 6) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
138 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
139 if (p.adj != "none") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
140 pvalue <- round(p.adjust(pvalue, p.adj), 6) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
141 LCL <- dif - Tprob * sdtdif |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
142 UCL <- dif + Tprob * sdtdif |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
143 sig <- rep(" ", nn) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
144 for (k in 1:nn) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
145 if (pvalue[k] <= 0.001) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
146 sig[k] <- "***" |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
147 else if (pvalue[k] <= 0.01) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
148 sig[k] <- "**" |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
149 else if (pvalue[k] <= 0.05) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
150 sig[k] <- "*" |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
151 else if (pvalue[k] <= 0.1) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
152 sig[k] <- "." |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
153 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
154 tr.i <- means[comb[1, ], 1] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
155 tr.j <- means[comb[2, ], 1] |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
156 dfComparisons <- data.frame(Difference = dif, p.value = pvalue, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
157 sig, LCL, UCL) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
158 rownames(dfComparisons) <- paste(tr.i, tr.j, sep = " - ") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
159 # cat("\nComparison between treatments mean of the ranks\n\n") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
160 # print(output) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
161 dfMeans <- data.frame(trt = means[, 1], means = means[, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
162 2], M = "", N = means[, 3]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
163 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
164 # invisible(output) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
165 invisible(list(study=main,test="Kruskal-Wallis test",value=H,df=(ntr - 1),chisq.p.value=p.chisq,p.adj.method=p.adj,ntStudent=dntStudent,alpha=alpha,LSD=dLSD,Harmonic.mean=dHMean,comparisons=dfComparisons,means=dfMeans)) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
166 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
167 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
168 ### This function is NOT original code but is from the gamlss package. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
169 ### It is written here in an effort to over write the gamlss object summary method |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
170 ### so that I can return information of interest. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
171 estimatesgamlss<-function (object, Qr, p1, coef.p, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
172 est.disp , df.r, |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
173 digits = max(3, getOption("digits") - 3), |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
174 covmat.unscaled , ...) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
175 { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
176 #covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE]) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
177 dimnames(covmat.unscaled) <- list(names(coef.p), names(coef.p)) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
178 covmat <- covmat.unscaled #in glm is=dispersion * covmat.unscaled, but here is already multiplied by the dispersion |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
179 var.cf <- diag(covmat) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
180 s.err <- sqrt(var.cf) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
181 tvalue <- coef.p/s.err |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
182 dn <- c("Estimate", "Std. Error") |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
183 if (!est.disp) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
184 { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
185 pvalue <- 2 * pnorm(-abs(tvalue)) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
186 coef.table <- cbind(coef.p, s.err, tvalue, pvalue) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
187 dimnames(coef.table) <- list(names(coef.p), c(dn, "z value","Pr(>|z|)")) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
188 } else if (df.r > 0) { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
189 pvalue <- 2 * pt(-abs(tvalue), df.r) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
190 coef.table <- cbind(coef.p, s.err, tvalue, pvalue) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
191 dimnames(coef.table) <- list(names(coef.p), c(dn, "t value","Pr(>|t|)")) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
192 } else { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
193 coef.table <- cbind(coef.p, Inf) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
194 dimnames(coef.table) <- list(names(coef.p), dn) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
195 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
196 return(coef.table) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
197 } |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
198 |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
199 ### This function is NOT original code but is from the gamlss package. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
200 ### It is written here in an effort to over write the gamlss object summary method |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
201 ### so that I can return information of interest. |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
202 summary.gamlss<- function (object, type = c("vcov", "qr"), save = FALSE, ...) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
203 { |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
204 return(as.data.frame(estimatesgamlss(object=object,Qr=object$mu.qr, p1=1:(object$mu.df-object$mu.nl.df), |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
205 coef.p=object$mu.coefficients[object$mu.qr$pivot[1:(object$mu.df-object$mu.nl.df)]], |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
206 est.disp =TRUE, df.r=(object$noObs - object$mu.df), |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
207 covmat.unscaled=chol2inv(object$mu.qr$qr[1:(object$mu.df-object$mu.nl.df), 1:(object$mu.df-object$mu.nl.df), drop = FALSE]) )) ) |
a87d5a5f2776
Uploaded the version running on the prod server
george-weingart
parents:
diff
changeset
|
208 } |