Guest User

Untitled

a guest
Jun 19th, 2018
91
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 8.23 KB | None | 0 0
  1. library(Hmisc)
  2. catTestchisq_sim_p <- function (tab)
  3. {
  4. st <- if (!is.matrix(tab) || nrow(tab) < 2 || ncol(tab) <
  5. 2)
  6. list(p.value = NA, statistic = NA, parameter = NA)
  7. else {
  8. rowcounts <- tab %*% rep(1, ncol(tab))
  9. tab <- tab[rowcounts > 0, ]
  10. if (!is.matrix(tab))
  11. list(p.value = NA, statistic = NA, parameter = NA)
  12. else chisq.test(tab, correct = FALSE, simulate.p.value = TRUE)
  13. }
  14. list(P = st$p.value, stat = st$statistic, df = st$parameter,
  15. testname = "Pearson", statname = "Chi-square", namefun = "chisq",
  16. latexstat = "\chi^{2}", plotmathstat = "chi^2")
  17. }
  18.  
  19. summaryM2 <- function (formula, groups = NULL, data = NULL, subset, na.action = na.retain,
  20. overall = FALSE, continuous = 10, na.include = FALSE, quant = c(0.025,
  21. 0.05, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 0.95,
  22. 0.975), nmin = 100, test = FALSE, conTest = conTestkw,
  23. catTest = catTestchisq_sim_p, ordTest = ordTestpo)
  24. {
  25. marg <- length(data) && ".marginal." %in% names(data)
  26. if (marg)
  27. formula <- update(formula, . ~ . + .marginal.)
  28. formula <- Formula(formula)
  29. Y <- if (!missing(subset) && length(subset))
  30. model.frame(formula, data = data, subset = subset, na.action = na.action)
  31. else model.frame(formula, data = data, na.action = na.action)
  32. X <- model.part(formula, data = Y, rhs = 1)
  33. Y <- model.part(formula, data = Y, lhs = 1)
  34. getlab <- function(x, default) {
  35. lab <- attr(x, "label")
  36. if (!length(lab) || lab == "")
  37. default
  38. else lab
  39. }
  40. if (marg) {
  41. xm <- X$.marginal.
  42. X$.marginal. <- NULL
  43. }
  44. else xm <- rep("", nrow(X))
  45. if (length(X)) {
  46. xname <- names(X)
  47. if (length(xname) == 1 && !length(groups))
  48. groups <- xname
  49. if (!length(groups) && length(xname) > 1) {
  50. warnings("Must specify groups when > 1 right hand side variable is present.ngroups taken as first right hand variable.")
  51. groups <- xname[1]
  52. }
  53. svar <- if (length(xname) == 1)
  54. factor(rep(".ALL.", nrow(X)))
  55. else do.call("interaction", list(X[setdiff(xname, groups)],
  56. sep = " "))
  57. group <- X[[groups]]
  58. glabel <- getlab(group, groups)
  59. }
  60. else {
  61. svar <- factor(rep(".ALL.", nrow(Y)))
  62. group <- rep("", nrow(Y))
  63. groups <- group.freq <- NULL
  64. glabel <- ""
  65. }
  66. quants <- unique(c(quant, 0.025, 0.05, 0.125, 0.25, 0.375,
  67. 0.5, 0.625, 0.75, 0.875, 0.95, 0.975))
  68. nv <- ncol(Y)
  69. nameY <- names(Y)
  70. R <- list()
  71. for (strat in levels(svar)) {
  72. instrat <- svar == strat
  73. n <- integer(nv)
  74. type <- n
  75. comp <- dat <- vector("list", nv)
  76. names(comp) <- names(dat) <- nameY
  77. labels <- Units <- vector("character", nv)
  78. if (test) {
  79. testresults <- vector("list", nv)
  80. names(testresults) <- names(comp)
  81. }
  82. gr <- group[instrat]
  83. xms <- xm[instrat]
  84. if (all(xms != ""))
  85. xms <- rep("", length(xms))
  86. group.freq <- table(gr)
  87. group.freq <- group.freq[group.freq > 0]
  88. if (overall)
  89. group.freq <- c(group.freq, Combined = sum(group.freq))
  90. for (i in 1:nv) {
  91. w <- Y[instrat, i]
  92. if (length(attr(w, "label")))
  93. labels[i] <- attr(w, "label")
  94. if (length(attr(w, "units")))
  95. Units[i] <- attr(w, "units")
  96. if (is.character(w))
  97. w <- as.factor(w)
  98. if (!inherits(w, "mChoice")) {
  99. if (!is.factor(w) && !is.logical(w) && length(unique(w[!is.na(w)])) <
  100. continuous)
  101. w <- as.factor(w)
  102. s <- !is.na(w)
  103. if (na.include && !all(s) && length(levels(w))) {
  104. w <- na.include(w)
  105. levels(w)[is.na(levels(w))] <- "NA"
  106. s <- rep(TRUE, length(s))
  107. }
  108. n[i] <- sum(s & xms == "")
  109. w <- w[s]
  110. g <- gr[s, drop = TRUE]
  111. if (is.factor(w) || is.logical(w)) {
  112. tab <- table(w, g)
  113. if (test) {
  114. if (is.ordered(w))
  115. testresults[[i]] <- ordTest(g, w)
  116. else testresults[[i]] <- catTest(tab)
  117. }
  118. if (nrow(tab) == 1) {
  119. b <- casefold(dimnames(tab)[[1]], upper = TRUE)
  120. pres <- c("1", "Y", "YES", "PRESENT")
  121. abse <- c("0", "N", "NO", "ABSENT")
  122. jj <- match(b, pres, nomatch = 0)
  123. if (jj > 0)
  124. bc <- abse[jj]
  125. else {
  126. jj <- match(b, abse, nomatch = 0)
  127. if (jj > 0)
  128. bc <- pres[jj]
  129. }
  130. if (jj) {
  131. tab <- rbind(tab, rep(0, ncol(tab)))
  132. dimnames(tab)[[1]][2] <- bc
  133. }
  134. }
  135. if (overall)
  136. tab <- cbind(tab, Combined = apply(tab, 1,
  137. sum))
  138. comp[[i]] <- tab
  139. type[i] <- 1
  140. }
  141. else {
  142. sfn <- function(x, quant) {
  143. y <- c(quantile(x, quant), Mean = mean(x),
  144. SD = sqrt(var(x)), N = sum(!is.na(x)))
  145. names(y) <- c(paste0(formatC(100 * quant,
  146. format = "fg", width = 1, digits = 15),
  147. "%"), "Mean", "SD", "N")
  148. y
  149. }
  150. qu <- tapply(w, g, sfn, simplify = TRUE, quants)
  151. if (test)
  152. testresults[[i]] <- conTest(g, w)
  153. if (overall)
  154. qu$Combined <- sfn(w, quants)
  155. comp[[i]] <- matrix(unlist(qu), ncol = length(quants) +
  156. 3, byrow = TRUE, dimnames = list(names(qu),
  157. c(format(quants), "Mean", "SD", "N")))
  158. if (any(group.freq <= nmin))
  159. dat[[i]] <- lapply(split(w, g), nmin = nmin,
  160. function(x, nmin) if (length(x) <= nmin)
  161. x
  162. else NULL)
  163. type[i] <- 2
  164. }
  165. }
  166. else {
  167. w <- as.numeric(w) == 1
  168. n[i] <- sum(!is.na(apply(w, 1, sum)) & xms ==
  169. "")
  170. g <- as.factor(gr)
  171. ncat <- ncol(w)
  172. tab <- matrix(NA, nrow = ncat, ncol = length(levels(g)),
  173. dimnames = list(dimnames(w)[[2]], levels(g)))
  174. if (test) {
  175. pval <- numeric(ncat)
  176. names(pval) <- dimnames(w)[[2]]
  177. d.f. <- stat <- pval
  178. }
  179. for (j in 1:ncat) {
  180. tab[j, ] <- tapply(w[, j], g, sum, simplify = TRUE,
  181. na.rm = TRUE)
  182. if (test) {
  183. tabj <- rbind(table(g) - tab[j, ], tab[j,
  184. ])
  185. st <- catTest(tabj)
  186. pval[j] <- st$P
  187. stat[j] <- st$stat
  188. d.f.[j] <- st$df
  189. }
  190. }
  191. if (test)
  192. testresults[[i]] <- list(P = pval, stat = stat,
  193. df = d.f., testname = st$testname, namefun = st$namefun,
  194. statname = st$statname, latexstat = st$latexstat,
  195. plotmathstat = st$plotmathstat)
  196. if (overall)
  197. tab <- cbind(tab, Combined = apply(tab, 1,
  198. sum))
  199. comp[[i]] <- tab
  200. type[i] <- 3
  201. }
  202. }
  203. labels <- ifelse(nchar(labels), labels, names(comp))
  204. R[[strat]] <- list(stats = comp, type = type, group.freq = group.freq,
  205. labels = labels, units = Units, quant = quant, data = dat,
  206. N = sum(!is.na(gr) & xms == ""), n = n, testresults = if (test) testresults)
  207. }
  208. structure(list(results = R, group.name = groups, group.label = glabel,
  209. call = call, formula = formula), class = "summaryM")
  210. }
  211.  
  212. options(digits=3)
  213. set.seed(173)
  214. sex <- factor(sample(c("m","f"), 500, rep=TRUE))
  215. treatment <- factor(sample(c("Drug","Placebo"), 500, rep=TRUE))
  216. f = summaryM2(sex ~ treatment, test=TRUE, overall = TRUE)
  217. latex(f, vnames = "labels",long = TRUE,
  218. exclude1 = FALSE, prmsd = FALSE, file = "", where = "!h", size = "scriptsize", digits = 2,
  219. what = "%", landscape = FALSE, longtable = TRUE, booktabs = TRUE, prtest= c('P', 'stat'),
  220. lines.page = (.Machine$integer.max + 9))
  221.  
  222. latex(f, prtest = c("P","stat","name"))
Add Comment
Please, Sign In to add comment