Meta analysis and clustering with special GMCMs"
In GMCM: Fast Estimation of Gaussian Mixture Copula Models

set.seed(5828993)
opts_knit$set(self.contained = TRUE)

Initialization

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.

Load data

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)") }}

Initial parameters

The initial parameters, as chosen in the application, are set to

init_par <- {{cput(init_par)}}

Model fitting

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}}.

Meta analysis with unsupervised clustering

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.

Results

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)))

Session information

This report was generated using rmarkdown^2^ and knitr^3^ under the session given below. The report utilizes parameterized reports and knitr::spin.

sessionInfo()

References

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"))