annotate src/lib/Misc.R @ 8:e9677425c6c3 default tip

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