-
Notifications
You must be signed in to change notification settings - Fork 2
/
meta_inla.R
203 lines (192 loc) · 8.72 KB
/
meta_inla.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
#' Fitting a pairwise meta-analysis model using INLA.
#'
#' \code{meta_inla} fits a pairwise meta-analysis model using INLA
#'
#' The following model types are supported
#' \itemize{
#' \item \code{FE}, fixed-effect model
#' \item \code{RE}, random effects model
#' }
#'
#' @param datINLA An object of \code{create_INLA_dat_pair}
#' @param fixed.par A numerical vector specifying the parameter of the normal prior
#' density for mean treatment effect, first value is parameter for mean, second is for variance.
#' @param tau.prior A string specifying the prior density for the heterogeneity standard deviation,
#' options are 'uniform' for uniform prior and 'half-normal' for half-normal prior.
#' @param tau.par A numerical vector specifying the parameter of the prior
#' density for heterogenety stdev.
#' \itemize{
#' \item \code{var.par = c(u, l)}: \code{u} is lower bound and \code{l} is upper
#' bound when \code{var.prior = 'uniform'}
#' \item \code{var.par = c(m, v)}: \code{m} is mean and \code{v} is variance
#' when \code{var.prior = 'uniform'}
#' }
#' @param mreg Logical indicating whether covariate(s) should be incorporated to fit a
#' meta-regression model, default \code{FALSE}
#' @param type A string indicating the type of the model, options are "FE", "RE".
#' @param approach A string indicating the approach of the model, options are "summary-level", "arm-level"
#' @param verbose Logical indicating whether the program should run in a verbose model, default \code{FALSE}.
#' @param inla.strategy A string specfying the strategy to use for the approximations of INLA;
#' one of 'gaussian', 'simplified.laplace' (default) or 'laplace', see \code{?INLA::control.inla}.
#' @param improve.hyperpar.dz Step length in the standardized scale used in the construction of the grid, default 0.75,
#' see \code{INLA::inla.hyperpar}.
#' @param correct Logical Add correction for the Laplace approximation, default \code{FALSE},
#' see \code{INLA::inla.hyperpar}.
#' @param correct.factor Numerical Factor used in adjusting the correction factor if \code{correct=TRUE}, default 10,
#' see \code{INLA::inla.hyperpar}.
#' @return \code{meta_inla} returns a \code{meta_inla} object with components:
#'
#' @examples
#' data('TBdat')
#' ## Create the dataset suitable for INLA
#' TBdatINLA <- create_INLA_dat_pair(TBdat$TRT, TBdat$CON, TBdat$TRTTB, TBdat$CONTB)
#'
#' ## Fitting a random-effects model using arm-level approach
#' \dontrun{
#' if(requireNamespace('INLA', quietly = TRUE)){
#' require('INLA', quietly = TRUE)
#' fit.TB.RE.INLA <- meta_inla(TBdatINLA, type = 'RE', approach = 'arm-level',
#' tau.prior = 'uniform', tau.par = c(0, 5))
#' }
#' }
#'
#' @export
meta_inla <- function(datINLA, fixed.par = c(0, 1000),
tau.prior = "uniform",
tau.par = c(0, 5),
type = "FE",
approach = "arm-level",
mreg = FALSE, verbose = FALSE, inla.strategy = "simplified.laplace", improve.hyperpar.dz = 0.75,
correct = FALSE, correct.factor = 10)
{
if (requireNamespace("INLA", quietly = TRUE)) {
if (!(sum(search() == "package:INLA")) == 1) {
stop("INLA need to be loaded! \n
Please use the following command to load INLA,\n
library(INLA) \n")
}
################ check data
if (!is.data.frame(datINLA$data.arm)) {
stop("Data MUST be a data frame!!!")
}
if(type %in% c("FE", "RE") == FALSE){
stop("Function argument \"type\" must be equal to \"FE\" or \"RE\"!")
}
if(approach %in% c("summary-level", "arm-level") == FALSE){
stop("Function argument \"approach\" must be equal to \"summary-level\" or \"arm-level\"!")
}
if(mreg == TRUE && is.null(datINLA$data.cont$cov)){
stop("Function argument \"cov\" must not be equal to \"NULL\" !")
}
inla.form <- "Y ~ -1 + d"
if(mreg == TRUE){
inla.form <- paste(inla.form, " + cov", sep = "")
}
if(approach == "arm-level"){
inla.form <- paste(inla.form, " + mu ", sep = "")
}
if(tau.prior == "half-normal") {
prior.expr <- " + f(het, model=\"iid\", prior = \"logtnormal\", param = tau.par"
}
if(tau.prior == "uniform") {
# Function for Uniform distribution:
hyperunif.function <- function(x) {
ifelse(exp(x)^-0.5 < tau.par[2] & exp(x)^-0.5 > tau.par[1], logdens <- log(1/(tau.par[2] - tau.par[1])), logdens <- log(9.98012604599318e-322))
logdenst <- logdens + log(0.5 * exp(-x/2))
return(logdenst)
}
# Set up grid to evaluate the uniform prior:
lprec <- seq(from = -40, to = 40, len = 20000) ## CHANGE this LINE if INLA crashes!
# (extend grid by changing (from= , to=))
# Create table with prior values and lprec:
unif.prior.table <- paste(c("table:", cbind(lprec, sapply(lprec, FUN = hyperunif.function))), sep = "", collapse = " ")
prior.expr <- " + f(het, model=\"iid\", hyper = list(theta = list(prior = unif.prior.table))"
}
if(type == "RE"){
inla.form <- paste(inla.form, prior.expr , ")", sep = "")
}
if(approach == "summary-level"){
prec = datINLA$data.cont$prec
fit.inla <- INLA::inla(stats::as.formula(inla.form), data = datINLA$data.cont, family = "normal",
control.fixed = list(expand.factor.strategy = "inla", mean = fixed.par[1],
prec = 1/fixed.par[2]), scale = prec,
control.family = list(hyper = list(prec = list( fixed = TRUE, initial = 0))),
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, mlik = TRUE, config = TRUE),
control.inla = list(strategy = inla.strategy, correct = correct, correct.factor = correct.factor))
if (!fit.inla$ok) {
stop("Something wrong while running model with data! Please set verbose = TRUE to check!!!!")
}
# improve the estimates for hyperparameters
fit.inla <- INLA::inla.hyperpar(fit.inla, dz = improve.hyperpar.dz)
}
if(approach == "arm-level"){
Ntrials = datINLA$data.arm$sampleSize
fit.inla <- INLA::inla(stats::as.formula(inla.form), data = datINLA$data.arm, family = "binomial",
control.fixed = list(expand.factor.strategy = "inla", mean = fixed.par[1],
prec = 1/fixed.par[2]), Ntrials = Ntrials,
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE, mlik = TRUE, config = TRUE),
control.inla = list(strategy = inla.strategy, control.correct = list(correct, correct.factor)))
if (!fit.inla$ok) {
stop("Something wrong while running model with data! Please set verbose = TRUE to check!!!!")
}
# improve the estimates for hyperparameters
fit.inla <- INLA::inla.hyperpar(fit.inla, dz = improve.hyperpar.dz)
}
nu <- as.matrix(fit.inla$summary.fixed[1, c(1, 2, 3, 4, 5)])
if(type == "RE"){
tau.mean <- INLA::inla.emarginal(function(x) 1/sqrt(x), fit.inla$marginals.hyperpar[["Precision for het"]])
tau <- INLA::inla.tmarginal(function(x) 1/sqrt(x), fit.inla$marginals.hyperpar[["Precision for het"]])
tau.stdev = sqrt(INLA::inla.emarginal(function(x) x^2, tau) - INLA::inla.emarginal(function(x) x^1, tau)^2)
tau.quant <- as.numeric(rev(sqrt((1/summary(fit.inla)$hyperpar[1, c(3, 4, 5)]))))
tab <- matrix(NA, 1, 5)
tab[1, ] <- c(tau.mean, tau.stdev, tau.quant[1], tau.quant[2], tau.quant[3])
colnames(tab) <- c("mean", "sd", "0.025quant", "0.5quant", "0.975quant")
rownames(tab) <- "tau"
} else tab <- NA
if(mreg == TRUE){
cov <- as.numeric(fit.inla$summary.fixed[2, c(1, 2, 3, 4, 5)])
}
fit.inla$nu <- nu
if (mreg == TRUE) {
fit.inla$cov <- cov
}
fit.inla$tau <- tab
fit.inla$fixed.par <- fixed.par
fit.inla$tau.prior <- tau.prior
fit.inla$tau.par <- tau.par
fit.inla$type <- type
fit.inla$approach <- approach
fit.inla$mreg <- mreg
fit.inla$inla.strategy <- inla.strategy
class(fit.inla) <- "meta_inla"
return(fit.inla)
} else {
stop("INLA need to be installed and loaded!\n
Please use the following command to install and load INLA,\n
install.packages(\"INLA\", repos=\"http://www.math.ntnu.no/inla/R/testing\") \n
library(INLA) \n")
}
}
#' @export
print.meta_inla <- function (x, digits = 2, ...)
{
if (!is.element("meta_inla", class(x)))
stop("Argument 'x' must be an object of class \"meta.inla\".")
cat("Time used: \n")
print(x$cpu.used)
if(x$mreg == FALSE){
cat("Meta analysis using INLA\n")
cat("Treatment effect\n")
print(round(x$nu, digits))
} else {
cat("Meta regression using INLA\n")
cat("Intercept\n")
print(round(x$nu, digits))
cat("Covariate\n")
print(round(x$cov, digits))
}
if(x$type == "RE"){
cat("Heterogeneity stdev\n")
print(round(x$tau, digits))
}
}