Estimate GMCM parameters of the special model

This function estimates the parameters of the special restricted Gaussian mixture copula model (GMCM) proposed by Li et. al. (2011). It is used to perform reproducibility (or meta) analysis using GMCMs. It features various optimization routines to identify the maximum likelihood estimate of the special GMCMs.

fit.meta.GMCM(u, init.par, method = c("NM", "SANN", "L-BFGS", "L-BFGS-B",
  "PEM"), max.ite = 1000, verbose = TRUE, positive.rho = TRUE,
  trace.theta = FALSE, ...)

fit.special.GMCM(u, init.par, method = c("NM", "SANN", "L-BFGS",
  "L-BFGS-B", "PEM"), max.ite = 1000, verbose = TRUE,
  positive.rho = TRUE, trace.theta = FALSE, ...)

Arguments

u	An `n` by `d` matrix of test statistics. Rows correspond to features and columns to experiments. Larger values are assumed to be indicative of stronger evidence and reproducibility.
init.par	A 4-dimensional vector of the initial parameters where, `init.par[1]` is the mixture proportion of spurious signals, `init.par[2]` is the mean, `init.par[3]` is the standard deviation, `init.par[4]` is the correlation.
method	A character vector of length \(1\). The optimization method used. Should be either `"NM"`, `"SANN"`, `"L-BFGS"`, `"L-BFGS-B"`, or `"PEM"` which are abbreviations of Nelder-Mead, Simulated Annealing, limited-memory quasi-Newton method, limited-memory quasi-Newton method with box constraints, and the pseudo EM algorithm, respectively. Default is `"NM"`. See `optim` for further details.
max.ite	The maximum number of iterations. If the `method` is `"SANN"` this is the number of iterations as there is no other stopping criterion. (See `optim`)
verbose	Logical. If `TRUE`, the log-likelihood values are printed.
positive.rho	`logical`. If `TRUE`, the correlation parameter is restricted to be positive.
trace.theta	`logical`. Extra convergence information is appended as a list to the output returned if `TRUE`. The exact behavior is dependent on the value of `method`. If `method` equals `"PEM"`, the argument is passed to `trace.theta` in `PseudoEMAlgorithm`. Otherwise it is passed to the control argument `trace` in `optim`.
...	Arguments passed to the `control`-list in `optim` or `PseudoEMAlgorithm` if `method` is `"PEM"`.

Value

A vector par of length 4 of the fitted parameters where par[1] is the probability of being from the first (or null) component, par[2] is the mean, par[3] is the standard deviation, and par[4] is the correlation.

If trace.theta is TRUE, then a list is returned where the first entry is as described above and the second entry is the trace information (dependent of method.).

Details

The "L-BFGS-B" method does not perform a transformation of the parameters.

fit.special.GMCM is simply an alias of fit.meta.gmcm.

Note

Simulated annealing is strongly dependent on the initial values and the cooling scheme.

See optim for further details.

References

Li, Q., Brown, J. B. J. B., Huang, H., & Bickel, P. J. (2011). Measuring reproducibility of high-throughput experiments. The Annals of Applied Statistics, 5(3), 1752-1779. doi:10.1214/11-AOAS466

Author

Anders Ellern Bilgrau <anders.ellern.bilgrau@gmail.com>

Examples

set.seed(1)

# True parameters
true.par <- c(0.9, 2, 0.7, 0.6)
# Simulation of data from the GMCM model
data <- SimulateGMCMData(n = 1000, par = true.par)
uhat <- Uhat(data$u) # Ranked observed data

init.par <- c(0.5, 1, 0.5, 0.9)  # Initial parameters

# Optimization with Nelder-Mead
nm.par   <- fit.meta.GMCM(uhat, init.par = init.par, method = "NM")
#>   Nelder-Mead direct search function minimizer
#> function value for initial parameters = 176.343401
#>   Scaled convergence tolerance is 2.62772e-06
#> Stepsize computed as 0.219722
#> BUILD              5 226.656977 105.973352
#> LO-REDUCTION       7 215.745222 105.973352
#> EXTENSION          9 176.343401 4.976728
#> LO-REDUCTION      11 133.560182 4.976728
#> EXTENSION         13 114.710466 -45.512205
#> EXTENSION         15 105.973352 -86.690136
#> EXTENSION         17 11.976424 -97.922305
#> LO-REDUCTION      19 4.976728 -97.922305
#> LO-REDUCTION      21 -45.512205 -97.922305
#> LO-REDUCTION      23 -78.907712 -97.922305
#> REFLECTION        25 -86.690136 -105.072127
#> LO-REDUCTION      27 -95.841590 -105.072127
#> LO-REDUCTION      29 -97.114144 -105.072127
#> LO-REDUCTION      31 -97.922305 -105.072127
#> LO-REDUCTION      33 -99.021762 -105.072127
#> LO-REDUCTION      35 -102.447782 -105.072127
#> EXTENSION         37 -102.566757 -109.005176
#> HI-REDUCTION      39 -103.544631 -109.005176
#> LO-REDUCTION      41 -105.010085 -109.005176
#> EXTENSION         43 -105.072127 -111.555517
#> LO-REDUCTION      45 -105.494019 -111.555517
#> EXTENSION         47 -107.926302 -114.706877
#> LO-REDUCTION      49 -108.163783 -114.706877
#> LO-REDUCTION      51 -109.005176 -114.706877
#> REFLECTION        53 -111.555517 -115.585400
#> HI-REDUCTION      55 -113.529332 -115.585400
#> LO-REDUCTION      57 -113.796380 -115.585400
#> EXTENSION         59 -113.979420 -116.852099
#> HI-REDUCTION      61 -114.706877 -116.852099
#> LO-REDUCTION      63 -115.194526 -116.852099
#> EXTENSION         65 -115.255631 -117.583240
#> HI-REDUCTION      67 -115.585400 -117.583240
#> EXTENSION         69 -116.276863 -118.192081
#> LO-REDUCTION      71 -116.284681 -118.192081
#> REFLECTION        73 -116.852099 -118.328874
#> LO-REDUCTION      75 -117.583240 -118.387443
#> LO-REDUCTION      77 -117.737665 -118.387443
#> HI-REDUCTION      79 -118.192081 -118.387443
#> HI-REDUCTION      81 -118.204023 -118.387443
#> HI-REDUCTION      83 -118.294593 -118.387443
#> LO-REDUCTION      85 -118.325318 -118.387443
#> REFLECTION        87 -118.328874 -118.401466
#> LO-REDUCTION      89 -118.331957 -118.401466
#> LO-REDUCTION      91 -118.368285 -118.401466
#> REFLECTION        93 -118.387443 -118.405678
#> EXTENSION         95 -118.396046 -118.435327
#> HI-REDUCTION      97 -118.398167 -118.435327
#> HI-REDUCTION      99 -118.401466 -118.435327
#> EXTENSION        101 -118.405678 -118.467822
#> LO-REDUCTION     103 -118.417378 -118.467822
#> LO-REDUCTION     105 -118.418910 -118.467822
#> EXTENSION        107 -118.435327 -118.480598
#> LO-REDUCTION     109 -118.440627 -118.480598
#> EXTENSION        111 -118.443666 -118.517242
#> EXTENSION        113 -118.467822 -118.589741
#> LO-REDUCTION     115 -118.478467 -118.589741
#> LO-REDUCTION     117 -118.480598 -118.589741
#> EXTENSION        119 -118.517242 -118.728493
#> LO-REDUCTION     121 -118.528818 -118.728493
#> EXTENSION        123 -118.534609 -118.882758
#> EXTENSION        125 -118.589741 -119.285949
#> LO-REDUCTION     127 -118.661470 -119.285949
#> EXTENSION        129 -118.728493 -119.665201
#> REFLECTION       131 -118.882758 -119.717315
#> LO-REDUCTION     133 -119.273816 -119.717315
#> LO-REDUCTION     135 -119.285949 -119.717315
#> HI-REDUCTION     137 -119.586776 -119.717315
#> HI-REDUCTION     139 -119.642925 -119.717315
#> EXTENSION        141 -119.662587 -119.756290
#> LO-REDUCTION     143 -119.665201 -119.756290
#> EXTENSION        145 -119.700026 -119.882576
#> LO-REDUCTION     147 -119.717315 -119.882576
#> EXTENSION        149 -119.751764 -120.001647
#> LO-REDUCTION     151 -119.756290 -120.001647
#> EXTENSION        153 -119.845613 -120.217221
#> LO-REDUCTION     155 -119.882576 -120.217221
#> LO-REDUCTION     157 -119.924760 -120.217221
#> EXTENSION        159 -120.001647 -120.535915
#> LO-REDUCTION     161 -120.176396 -120.535915
#> EXTENSION        163 -120.209100 -120.720566
#> LO-REDUCTION     165 -120.217221 -120.720566
#> EXTENSION        167 -120.525797 -121.072510
#> LO-REDUCTION     169 -120.535915 -121.072510
#> LO-REDUCTION     171 -120.710482 -121.072510
#> REFLECTION       173 -120.720566 -121.132790
#> HI-REDUCTION     175 -120.974799 -121.132790
#> HI-REDUCTION     177 -120.983393 -121.132790
#> LO-REDUCTION     179 -121.034027 -121.132790
#> LO-REDUCTION     181 -121.059421 -121.132790
#> HI-REDUCTION     183 -121.072510 -121.132790
#> LO-REDUCTION     185 -121.089634 -121.132790
#> EXTENSION        187 -121.103641 -121.188826
#> LO-REDUCTION     189 -121.111154 -121.188826
#> LO-REDUCTION     191 -121.122667 -121.188826
#> EXTENSION        193 -121.132790 -121.204569
#> EXTENSION        195 -121.134706 -121.222492
#> EXTENSION        197 -121.180724 -121.336842
#> LO-REDUCTION     199 -121.188826 -121.336842
#> LO-REDUCTION     201 -121.204569 -121.336842
#> LO-REDUCTION     203 -121.222492 -121.336842
#> REFLECTION       205 -121.282353 -121.360958
#> EXTENSION        207 -121.283241 -121.432059
#> EXTENSION        209 -121.329228 -121.555775
#> LO-REDUCTION     211 -121.336842 -121.555775
#> LO-REDUCTION     213 -121.360958 -121.555775
#> EXTENSION        215 -121.432059 -121.725111
#> LO-REDUCTION     217 -121.444713 -121.725111
#> EXTENSION        219 -121.529250 -121.904929
#> EXTENSION        221 -121.555775 -122.133960
#> EXTENSION        223 -121.692835 -122.326879
#> LO-REDUCTION     225 -121.725111 -122.326879
#> LO-REDUCTION     227 -121.904929 -122.326879
#> REFLECTION       229 -122.133960 -122.383178
#> HI-REDUCTION     231 -122.199466 -122.383178
#> EXTENSION        233 -122.211678 -122.561704
#> LO-REDUCTION     235 -122.282918 -122.561704
#> LO-REDUCTION     237 -122.326879 -122.561704
#> REFLECTION       239 -122.383178 -122.648447
#> REFLECTION       241 -122.541034 -122.732974
#> HI-REDUCTION     243 -122.557909 -122.732974
#> LO-REDUCTION     245 -122.561704 -122.732974
#> REFLECTION       247 -122.626764 -122.761473
#> EXTENSION        249 -122.648447 -122.806684
#> LO-REDUCTION     251 -122.689278 -122.806684
#> REFLECTION       253 -122.732974 -122.824853
#> LO-REDUCTION     255 -122.761473 -122.824853
#> LO-REDUCTION     257 -122.782913 -122.824877
#> HI-REDUCTION     259 -122.806684 -122.824877
#> LO-REDUCTION     261 -122.809614 -122.824877
#> LO-REDUCTION     263 -122.824241 -122.830195
#> HI-REDUCTION     265 -122.824791 -122.830195
#> LO-REDUCTION     267 -122.824853 -122.830891
#> HI-REDUCTION     269 -122.824877 -122.830891
#> LO-REDUCTION     271 -122.829721 -122.830891
#> HI-REDUCTION     273 -122.830143 -122.831148
#> LO-REDUCTION     275 -122.830195 -122.831543
#> HI-REDUCTION     277 -122.830342 -122.831543
#> HI-REDUCTION     279 -122.830891 -122.831543
#> HI-REDUCTION     281 -122.831148 -122.831742
#> LO-REDUCTION     283 -122.831400 -122.831760
#> REFLECTION       285 -122.831449 -122.831854
#> HI-REDUCTION     287 -122.831543 -122.831854
#> REFLECTION       289 -122.831742 -122.831879
#> LO-REDUCTION     291 -122.831760 -122.831916
#> HI-REDUCTION     293 -122.831785 -122.831918
#> REFLECTION       295 -122.831854 -122.832122
#> HI-REDUCTION     297 -122.831879 -122.832122
#> LO-REDUCTION     299 -122.831916 -122.832122
#> HI-REDUCTION     301 -122.831918 -122.832134
#> LO-REDUCTION     303 -122.832057 -122.832134
#> LO-REDUCTION     305 -122.832059 -122.832134
#> EXTENSION        307 -122.832103 -122.832258
#> LO-REDUCTION     309 -122.832117 -122.832258
#> LO-REDUCTION     311 -122.832122 -122.832258
#> HI-REDUCTION     313 -122.832134 -122.832258
#> LO-REDUCTION     315 -122.832189 -122.832258
#> EXTENSION        317 -122.832225 -122.832323
#> LO-REDUCTION     319 -122.832237 -122.832323
#> LO-REDUCTION     321 -122.832258 -122.832323
#> REFLECTION       323 -122.832258 -122.832331
#> LO-REDUCTION     325 -122.832304 -122.832331
#> LO-REDUCTION     327 -122.832308 -122.832331
#> LO-REDUCTION     329 -122.832320 -122.832334
#> REFLECTION       331 -122.832323 -122.832336
#> EXTENSION        333 -122.832325 -122.832341
#> REFLECTION       335 -122.832331 -122.832343
#> REFLECTION       337 -122.832334 -122.832346
#> EXTENSION        339 -122.832336 -122.832357
#> LO-REDUCTION     341 -122.832341 -122.832357
#> REFLECTION       343 -122.832343 -122.832361
#> LO-REDUCTION     345 -122.832346 -122.832361
#> REFLECTION       347 -122.832352 -122.832362
#> LO-REDUCTION     349 -122.832354 -122.832362
#> REFLECTION       351 -122.832357 -122.832364
#> REFLECTION       353 -122.832361 -122.832365
#> HI-REDUCTION     355 -122.832362 -122.832365
#> HI-REDUCTION     357 -122.832362 -122.832365
#> Exiting from Nelder Mead minimizer
#>     359 function evaluations used

if (FALSE) {
# Comparison with other optimization methods
# Optimization with simulated annealing
sann.par <- fit.meta.GMCM(uhat, init.par = init.par, method = "SANN",
                          max.ite = 3000, temp = 1)
# Optimization with the Pseudo EM algorithm
pem.par  <- fit.meta.GMCM(uhat, init.par = init.par, method = "PEM")

# The estimates agree nicely
rbind("True" = true.par, "Start" = init.par,
      "NM" = nm.par, "SANN" = sann.par, "PEM" = pem.par)
}

# Get estimated cluster
Khat <- get.IDR(x = uhat, par = nm.par)$Khat
plot(uhat, col = Khat, main = "Clustering\nIDR < 0.05")