set.seed(5828993) opts_knit$set(self.contained = TRUE)
The GMCM^1^ package is loaded.
#install.packages("GMCM") # Uncomment to install the GMCM package library("GMCM")
If GMCM is not installed, please uncomment the above line and rerun the script.
The data is loaded and the first rows are printed
ds <- read.table(file = "{{gsub("\\", "/", data_file, fixed = TRUE)}}", header = {{header}}, sep = "{{sep}}", quote = "\{{quote}}") head(ds, n = 4)
To reproduce the output, change the file path above to the data file of interest.
Next, the data is subset to the columns of interest.
x <- ds[, {{cput(model_cols)}}] head(x, n = 2)
The rescaling is modified if the 'large values indicate strong evidence' check box was checked or not.
The check box was {{ ifelse(meta_large_vals, "checked.", "unchecked, so we reverse ranking.") }}
{{ ifelse(meta_large_vals, "u <- Uhat(x)", "u <- Uhat(-x)") }}
The initial parameters, as chosen in the application, are set to
init_par <- {{cput(init_par)}}
With the data loaded and defined initial parameters, the model is can be fitted.
par <- fit.meta.GMCM(u = u, init.par = init_par, method = "{{meta_method}}", max.ite = {{meta_max_ite}}, positive.rho = {{meta_positive_rho}}, verbose = FALSE) print(par)
I.e. we set the fitting method to "{{meta_method}}"
with a maximum number of iterations of {{meta_max_ite}}
.
The estimated parameters are used to calculate the local and adjusted irreproducibility discovery rates:
meta_IDR_thres <- {{meta_IDR_thres}} out <- get.IDR(x, par = par, threshold = meta_IDR_thres) # Compute IDR str(out) # Use idr instead of IDR to threshold if specified out <- get.IDR(u, par = par, threshold = meta_IDR_thres) below <- (out[["{{meta_IDR_thres_type}}"]] < meta_IDR_thres) out$l <- sum(below) out$Khat <- ifelse(below, 2, 1)
By default, get.IDR
computes thresholds on the adjusted IDR.
The local irreproducibility discovery rate correspond to the posterior probability of the point originating from the irreproducible component.
The classes are counted by
table(out$Khat)
The results are also displayed by plotting
plot(x, col = out$Khat, asp = 1) # Plot of raw values plot(u, col = out$Khat, asp = 1) # Plot of copula values z <- GMCM:::qgmm.marginal(u, theta = meta2full(par, d = ncol(u))) # Estimate latent process plot(z, col = out$Khat, asp = 1) # Plot of estimated latent process
The fitted par
object can be converted to a theta
object which can be plotted directly:
plot(meta2full(par, d = ncol(u)))
This report was generated using rmarkdown^2^ and knitr^3^ under the session
given below. The report utilizes parameterized reports and knitr::spin
.
sessionInfo()
Please cite the GMCM paper^1^ if you use the package or shiny app.
cites <- lapply(c("GMCM", "knitr", "rmarkdown"), citation) fmt_cites <- unlist(lapply(cites, format, style = "text")) cat(paste0(" ", seq_along(fmt_cites), ". ", fmt_cites, "\n"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.