Function to profile the (restricted) log-likelihood for objects of class "rma.uni", "rma.mv", "rma.uni.selmodel", and "rma.ls".

# S3 method for rma.uni
profile(fitted, xlim, ylim, steps=20,
        lltol=1e-03, progbar=TRUE, parallel="no", ncpus=1, cl=NULL,
        plot=TRUE, pch=19, refline=TRUE, cline=FALSE, ...)

# S3 method for rma.mv
profile(fitted, sigma2, tau2, rho, gamma2, phi,
        xlim, ylim, steps=20, lltol=1e-03, progbar=TRUE,
        parallel="no", ncpus=1, cl=NULL,
        plot=TRUE, pch=19, refline=TRUE, cline=FALSE, ...)

# S3 method for rma.uni.selmodel
profile(fitted, tau2, delta,
        xlim, ylim, steps=20, lltol=1e-03, progbar=TRUE,
        parallel="no", ncpus=1, cl=NULL,
        plot=TRUE, pch=19, refline=TRUE, cline=FALSE, ...)

# S3 method for rma.ls
profile(fitted, alpha, xlim, ylim, steps=20,
        lltol=1e-03, progbar=TRUE, parallel="no", ncpus=1, cl=NULL,
        plot=TRUE, pch=19, refline=TRUE, cline=FALSE, ...)

# S3 method for profile.rma
print(x, ...)
# S3 method for profile.rma
plot(x, xlim, ylim, pch=19,
     xlab, ylab, main, refline=TRUE, cline=FALSE, ...)

Arguments

fitted

an object of class "rma.uni", "rma.mv", "rma.uni.selmodel", or "rma.ls".

x

an object of class "profile.rma" (for plot and print).

sigma2

optional integer to specify for which \(\sigma^2\) value the likelihood should be profiled.

tau2

optional integer to specify for which \(\tau^2\) value the likelihood should be profiled.

rho

optional integer to specify for which \(\rho\) value the likelihood should be profiled.

gamma2

optional integer to specify for which \(\gamma^2\) value the likelihood should be profiled.

phi

optional integer to specify for which \(\phi\) value the likelihood should be profiled.

delta

optional integer to specify for which \(\delta\) value the likelihood should be profiled.

alpha

optional integer to specify for which \(\alpha\) value the likelihood should be profiled.

xlim

optional vector to specify the lower and upper limit of the parameter over which the profiling should be done. If unspecified, the function tries to set these limits automatically.

ylim

optional vector to specify the y-axis limits when plotting the profiled likelihood. If unspecified, the function tries to set these limits automatically.

steps

number of points between xlim[1] and xlim[2] (inclusive) for which the likelihood should be evaluated (the default is 20).

lltol

numerical tolerance used when comparing values of the profiled log-likelihood with the log-likelihood of the fitted model (the default is 1e-03).

progbar

logical to specify whether a progress bar should be shown (the default is TRUE).

parallel

character string to specify whether parallel processing should be used (the default is "no"). For parallel processing, set to either "snow" or "multicore". See ‘Details’.

ncpus

integer to specify the number of processes to use in the parallel processing.

cl

optional cluster to use if parallel="snow". If not supplied, a cluster on the local machine is created for the duration of the call.

plot

logical to specify whether the profile plot should be drawn after profiling is finished (the default is TRUE).

pch

plotting symbol to use. By default, a filled circle is used. See points for other options.

refline

logical to specify whether the value of the parameter estimate should be indicated by a dotted vertical line and its log-likelihood value by a dotted horizontal line (the default is TRUE).

cline

logical to specify whether a horizontal reference line should be added to the plot that indicates the log-likelihood value corresponding to the 95% profile confidence interval (the default is FALSE).

xlab

title for the x-axis. If unspecified, the function tries to set an appropriate axis title.

ylab

title for the y-axis. If unspecified, the function tries to set an appropriate axis title.

main

title for the plot. If unspecified, the function tries to set an appropriate title.

...

other arguments.

Details

The function fixes a particular parameter of the model and then computes the maximized (restricted) log-likelihood over the remaining parameters of the model. By doing this for a range of values for the parameter that was fixed, a profile of the (restricted) log-likelihood is constructed.

The parameter(s) that can be profiled depend on the model object:

  • For objects of class "rma.uni" obtained with the rma.uni function, the function profiles over parameter \(\tau^2\) (not for fixed-effects models).

  • For objects of class "rma.mv" obtained with the rma.mv function, profiling is done by default over all (non-fixed) variance and correlation components of the model. Alternatively, one can use the sigma2, tau2, rho, gamma2, or phi arguments to specify over which parameter the profiling should be done. Only one of these arguments can be used at a time. A single integer is used to specify the number of the parameter.

  • For selection model objects of class "rma.uni.selmodel" obtained with the selmodel.rma.uni function, profiling is done by default over \(\tau^2\) (for models where this is an estimated parameter) and all (non-fixed) selection model parameters. Alternatively, one can choose to profile only \(\tau^2\) by setting tau2=TRUE or one can select one of the selection model parameters to profile by specifying its number via the delta argument.

  • For location-scale model objects of class "rma.ls" obtained with the rma.uni function, profiling is done by default over all (non-fixed) \(\alpha\) parameters that are part of the scale model.

A profile plot should show a single peak at the corresponding ML/REML estimate (assuming that the model was fitted with ML/REML estimation). If refline=TRUE (the default), the value of the parameter estimate is indicated by a dotted vertical line and its log-likelihood value by a dotted horizontal line. Hence, the intersection of these two lines should correspond to the peak.

When profiling a variance component (or some other parameter that cannot be negative), the peak may be at zero (if this corresponds to the ML/REML estimate of the parameter). In this case, the profiled log-likelihood should be a monotonically decreasing function of the parameter.

If the profiled log-likelihood has multiple peaks, this indicates that the likelihood surface is not unimodal. In such cases, the ML/REML estimate may correspond to a local optimum (when the intersection of the two dotted lines is not at the highest peak).

If the profile is flat (over the entire parameter space or large portions of it), then this suggests that at least some of the parameters of the model are not identifiable (and the parameter estimates obtained are to some extent arbitrary). See Raue et al. (2009) for some further discussion of parameter identifiability (structurally and practically) and the use of profile likelihoods to check for this.

The function checks whether any profiled log-likelihood value is actually larger than the log-likelihood of the fitted model (using a numerical tolerance of lltol). If so, a warning is issued as this might indicate that the optimizer did not identify the actual ML/REML estimate.

Profiling requires repeatedly refitting the model, which can be slow when \(k\) is large and/or the model is complex (the latter especially applies to "rma.mv" objects and also to certain "rma.uni.selmodel" or "rma.ls" objects). On machines with multiple cores, one can usually speed things up by delegating the model fitting to separate worker processes, that is, by setting parallel="snow" or parallel="multicore" and ncpus to some value larger than 1. Parallel processing makes use of the parallel package, using the makePSOCKcluster and parLapply functions when parallel="snow" or using mclapply when parallel="multicore" (the latter only works on Unix/Linux-alikes). With parallel::detectCores(), one can check on the number of available cores on the local machine.

Value

An object of class "profile.rma". The object is a list (or list of lists) containing the following components:

sigma2

values of \(\sigma^2\) over which the likelihood was profiled (only when profiling was actually done over \(\sigma^2\)).

tau2

values of \(\tau^2\) over which the likelihood was profiled (only when profiling was actually done over \(\tau^2\)).

rho

values of \(\rho\) over which the likelihood was profiled (only when profiling was actually done over \(\rho\)).

gamma2

values of \(\gamma^2\) over which the likelihood was profiled (only when profiling was actually done over \(\gamma^2\)).

phi

values of \(\phi\) over which the likelihood was profiled (only when profiling was actually done over \(\phi\)).

delta

values of \(\delta\) over which the likelihood was profiled (only when profiling was actually done over \(\delta\)).

alpha

values of \(\alpha\) over which the likelihood was profiled (only when profiling was actually done over \(\alpha\)).

ll

(restricted) log-likelihood at the corresponding parameter value.

beta

a matrix with the estimated model coefficients at the corresponding parameter value.

ci.lb

a matrix with the lower confidence interval bounds for the model coefficients at the corresponding parameter value.

ci.ub

a matrix with the upper confidence interval bounds for the model coefficients at the corresponding parameter value.

...

some additional elements/values.

Note that the list is returned invisibly.

Author

Wolfgang Viechtbauer wvb@metafor-project.org https://www.metafor-project.org

References

Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmuller, U., & Timmer, J. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923--1929. https://doi.org/10.1093/bioinformatics/btp358

Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1--48. https://doi.org/10.18637/jss.v036.i03

See also

Examples

### calculate log odds ratios and corresponding sampling variances dat <- escalc(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg) ### fit random-effects model using rma.uni() res <- rma(yi, vi, data=dat) ### profile over tau^2 profile(res, progbar=FALSE)
### change data into long format dat.long <- to.long(measure="OR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg) ### set levels of group variable ("exp" = experimental/vaccinated; "con" = control/non-vaccinated) levels(dat.long$group) <- c("exp", "con") ### set "con" to reference level dat.long$group <- relevel(dat.long$group, ref="con") ### calculate log odds and corresponding sampling variances dat.long <- escalc(measure="PLO", xi=out1, mi=out2, data=dat.long) ### fit bivariate random-effects model using rma.mv() res <- rma.mv(yi, vi, mods = ~ group, random = ~ group | study, struct="UN", data=dat.long) res
#> #> Multivariate Meta-Analysis Model (k = 26; method: REML) #> #> Variance Components: #> #> outer factor: study (nlvls = 13) #> inner factor: group (nlvls = 2) #> #> estim sqrt k.lvl fixed level #> tau^2.1 1.5486 1.2444 13 no con #> tau^2.2 2.6173 1.6178 13 no exp #> #> rho.con rho.exp con exp #> con 1 - 13 #> exp 0.9450 1 no - #> #> Test for Residual Heterogeneity: #> QE(df = 24) = 5270.3863, p-val < .0001 #> #> Test of Moderators (coefficient 2): #> QM(df = 1) = 15.5470, p-val < .0001 #> #> Model Results: #> #> estimate se zval pval ci.lb ci.ub #> intrcpt -4.8374 0.3528 -13.7113 <.0001 -5.5289 -4.1459 *** #> groupexp 0.7414 0.1880 3.9430 <.0001 0.3729 1.1099 *** #> #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #>
### profile over tau^2_1, tau^2_2, and rho ### note: for rho, adjust region over which profiling is done ('zoom in' on area around estimate) # \dontrun{ par(mfrow=c(3,1)) profile(res, tau2=1) profile(res, tau2=2) profile(res, rho=1, xlim=c(.90, .98))# }
### an example where the peak is at 0 dat <- escalc(measure="RD", n1i=n1i, n2i=n2i, ai=ai, ci=ci, data=dat.hine1989) res <- rma(yi, vi, data=dat) par(mfrow=c(1,1)) profile(res, progbar=FALSE)