1
#' @title geeUni: The coefficient of univariate gee (internal function)
2
#' @description Extract the coefficients of univariate gee using geeglm function (geepack package).
3
#' @param y Dependant variable
4
#' @param x Independent variable
5
#' @param data Data
6
#' @param id.vec Vector of id (should be ordered)
7
#' @param family Family: "gaussian", "binomial", "poisson", "quasipoisson", etc...
8
#' @param cor.type Correlation structure, Default: 'exchangeable'
9
#' @return coefficient, standard error, p-value
10
#' @details DETAILS
11
#' @examples 
12
#'  library(geepack)
13
#'  data(dietox)
14
#'  dietox$Cu <- as.factor(dietox$Cu)
15
#'  gee.uni <- geeUni("Weight", "Time", data = dietox, id.vec = dietox$Pig, 
16
#'                    family = "gaussian", cor.type = "exchangeable")
17
#' @rdname geeUni
18
#' @importFrom geepack geeglm 
19
#' @importFrom stats as.formula
20
#' @export
21

22

23
geeUni = function(y, x, data, id.vec, family, cor.type = "exchangeable"){
24 1
  form <- as.formula(paste(y, "~", x))
25 1
  res <- geepack::geeglm(form, data = data, family = family, id = id.vec, corstr = cor.type)
26 1
  coef <- summary(res)$coefficients[-1, -3]
27 1
  return(coef)
28
}
29

30

31

32
#' @title geeExp: transform the unit of coefficients (internal function)
33
#' @description Transform the unit of coefficients to "Coeff", "OR" or "RR"
34
#' @param gee.coef geeUni object.
35
#' @param family Family: "gaussian", "binomial", "poisson", "quasipoisson", etc..., Default: 'binomial'
36
#' @param dec Decimal point
37
#' @return The transforemed coefficients(95% CI), p-value
38
#' @details DETAILS
39
#' @examples 
40
#'  library(geepack)
41
#'  data(dietox)
42
#'  dietox$Cu <- as.factor(dietox$Cu)
43
#'  gee.uni <- geeUni("Weight", c("Time", "Cu"), data = dietox, id.vec = dietox$Pig, 
44
#'                    family = "gaussian", cor.type = "exchangeable")
45
#'  gee.exp <- geeExp(gee.uni, "binomial", 2)
46
#' @rdname geeExp
47
#' @export
48

49

50
geeExp = function(gee.coef, family ="binomial", dec){
51 1
  if (family == "binomial"){
52 1
    OR <- paste(round(exp(gee.coef[,1]), dec), " (", round(exp(gee.coef[,1] - 1.96*gee.coef[,2]), dec), ",", round(exp(gee.coef[,1] + 1.96*gee.coef[,2]), dec),")", sep="")
53 1
    return(cbind(OR, gee.coef[,3]))
54 1
  } else if (family == "gaussian"){
55 1
    coeff <- paste(round(gee.coef[,1], dec), " (", round(gee.coef[,1] - 1.96*gee.coef[,2], dec), ",", round(gee.coef[,1] + 1.96*gee.coef[,2], dec), ")", sep="")
56 1
    return(cbind(coeff, gee.coef[,3]))
57 0
  } else if (family %in% c("poisson", "quasipoisson")){
58 0
    RR <- paste(round(exp(gee.coef[,1]), dec), " (", round(exp(gee.coef[,1] - 1.96*gee.coef[,2]), dec), ",", round(exp(gee.coef[,1] + 1.96*gee.coef[,2]), dec),")", sep="")
59 0
    return(cbind(RR, gee.coef[,3]))
60
  }
61
} 
62

63

64

65
#' @title geeglm.display
66
#' @description Make gee results from "geeglm" object
67
#' @param geeglm.obj "geeglm" object
68
#' @param decimal Decimal, Default: 2
69
#' @return List: caption, main table, metrics table
70
#' @details DETAILS
71
#' @examples 
72
#'  library(geepack)
73
#'  data(dietox)
74
#'  dietox$Cu <- as.factor(dietox$Cu)
75
#'  gee01 <- geeglm (Weight ~ Time + Cu , id =Pig, data = dietox,
76
#'                 family=gaussian,corstr="ex")
77
#'  geeglm.display(gee01)
78
#' @seealso 
79
#'  \code{\link[data.table]{data.table-package}}
80
#'  \code{\link[stats]{complete.cases}}
81
#' @rdname geeglm.display
82
#' @export 
83
#' @importFrom data.table data.table
84
#' @importFrom stats complete.cases
85

86
geeglm.display = function(geeglm.obj, decimal = 2){
87 1
  family.gee <- geeglm.obj$family[[1]]
88 1
  corstr.gee <- geeglm.obj$corstr
89 1
  y <- as.character(geeglm.obj$terms[[2]])
90 1
  xs <- names(geeglm.obj$model)[-1]
91
  ## rownames
92 1
  geeglm.obj$data <- data.table::data.table(geeglm.obj$data)
93 1
  nomiss <- stats::complete.cases(geeglm.obj$data[, c(y, xs), with = F])
94 1
  gee.uni.list <- lapply(xs, function(x){geeUni(y, x, data = geeglm.obj$data[nomiss, ], id.vec = geeglm.obj$id, family = family.gee, cor.type = corstr.gee)})
95 1
  rn.uni <- lapply(gee.uni.list, rownames)
96 1
  gee.uni <- Reduce(rbind, gee.uni.list)
97
  
98 1
  if (length(xs) ==1){
99 1
    gee.res <- geeExp(gee.uni, family = family.gee, dec = decimal)
100 1
    gee.res.list <- lapply(1:length(xs), function(x){gee.res[rownames(gee.uni) %in% rn.uni[[x]], ]})     
101 1
    varnum.mfac <- which(lapply(gee.res.list, length) > ncol(gee.res))
102 0
    lapply(varnum.mfac, function(x){gee.res.list[[x]] <<- rbind(rep(NA, ncol(gee.res)), gee.res.list[[x]])})
103 1
    gee.res.modi <- Reduce(rbind, gee.res.list)
104 1
    if (nrow(gee.uni) == 1){
105 1
      gee.res.modi <- t(data.frame(gee.res.modi))
106
    } 
107

108 1
    family.label <- colnames(gee.res.modi)[1]
109 1
    colnames(gee.res.modi) <-  c(paste(family.label, "(95%CI)",sep = ""), "P value")
110
  } else{
111 1
    gee.multi <- summary(geeglm.obj)$coefficients[-1, -3]
112 1
    gee.res <- cbind(geeExp(gee.uni, family = family.gee, dec = decimal), geeExp(gee.multi, family = family.gee, dec = decimal))
113 1
    gee.res.list <- lapply(1:length(xs), function(x){gee.res[rownames(gee.uni) %in% rn.uni[[x]], ]})      
114 1
    varnum.mfac <- which(lapply(gee.res.list, length) > ncol(gee.res))
115 1
    lapply(varnum.mfac, function(x){gee.res.list[[x]] <<- rbind(rep(NA, ncol(gee.res)), gee.res.list[[x]])})
116 1
    gee.res.modi <- Reduce(rbind, gee.res.list)
117 1
    family.label <- colnames(gee.res.modi)[1]
118 1
    colnames(gee.res.modi) <- c(paste("crude ", family.label, "(95%CI)",sep = ""), "crude P value", paste("adj. ", family.label, "(95%CI)",sep = ""), "adj. P value")
119
  }
120
  
121
  
122
  
123 1
  rn.list <- lapply(1:length(xs), function(x){rownames(gee.uni)[rownames(gee.uni) %in% rn.uni[[x]]]})
124 1
  varnum.2fac <- which(lapply(xs, function(x){length(geeglm.obj$xlevels[[x]])}) == 2)
125 0
  lapply(varnum.2fac, function(x){rn.list[[x]] <<- paste(xs[x], ": ", geeglm.obj$xlevels[[xs[x]]][2], " vs ", geeglm.obj$xlevels[[xs[x]]][1], sep="")})
126 1
  lapply(varnum.mfac, function(x){rn.list[[x]] <<- c(paste(xs[x],": ref.=", geeglm.obj$xlevels[[xs[x]]][1], sep=""), gsub(xs[x],"   ", rn.list[[x]]))})
127
  
128
  
129 1
  rownames(gee.res.modi) <- unlist(rn.list)
130
  #out = as.data.frame(gee.res.modi)
131
  #lapply(seq(2, ncol(gee.res), by =2), function(x){out[, x] <<- as.numeric(as.vector(out[, x]))})
132 1
  out <- gee.res.modi
133 1
  pv.colnum <- which(colnames(out) %in% c("P value", "crude P value", "adj. P value"))
134 1
  for (i in pv.colnum){
135 1
    out[, i] <- ifelse(as.numeric(out[, i]) < 0.001, "< 0.001", round(as.numeric(out[, i]), decimal + 1))
136
  }
137
  
138
  ## Metric
139 1
  info.gee <- as.character(c(NA, round(as.numeric(summary(geeglm.obj)$corr[1]), decimal+1), length(unique(geeglm.obj$id)), length(geeglm.obj$y)))
140 1
  info.df <- cbind(info.gee, matrix(NA, 4, ncol(gee.res) -1))
141 1
  colnames(info.df) <- colnames(out)
142
  #lapply(seq(2, ncol(gee.res), by =2), function(x){info.df[, x] <<- as.numeric(info.df[, x])})
143 1
  rownames(info.df) <- c("","Estimated correlation parameters", "No. of clusters", "No. of observations")
144
  
145
  ## Caption
146 1
  cap.gee = paste("GEE(", family.gee, ") predicting ", y, " by ", paste(xs, collapse = ", "), " - Group ", as.character(geeglm.obj$call$id)[length(as.character(geeglm.obj$call$id))], sep="")
147
  
148 1
  return(list(caption = cap.gee, table = out, metric = info.df))
149
}

Read our documentation on viewing source code .

Loading