Writes MLwiN macros to fit models using Markov chain Monte Carlo (MCMC) methods

Share:

Description

write.MCMC is an internal function which creates an MLwiN macro file to fit models using MCMC.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
write.MCMC(indata, dtafile, oldsyntax = FALSE, resp, levID, expl, rp,
  D = "Normal", nonlinear = c(0, 1), categ = NULL, notation = NULL,
  nonfp = NULL, clre = NULL, Meth = 1, merr = NULL, carcentre = FALSE,
  maxiter = 20, convtol = 2, seed = 1, iterations = 5000,
  burnin = 500, scale = 5.8, thinning = 1, priorParam = "default",
  refresh = 50, fixM = 1, residM = 1, Lev1VarM = 1, OtherVarM = 1,
  adaption = 1, priorcode = c(gamma = 1), rate = 50, tol = 10,
  lclo = 0, mcmcOptions, fact = NULL, xc = NULL, mm = NULL,
  car = NULL, BUGO = NULL, mem.init = "default", optimat = FALSE,
  modelfile, initfile, datafile, macrofile, IGLSfile, MCMCfile, chainfile,
  MIfile, resifile, resi.store = FALSE, resioptions, resichains,
  FACTchainfile, resi.store.levs = NULL, debugmode = FALSE,
  startval = NULL, dami = NULL, namemap = sapply(colnames(indata),
  as.character), saveworksheet = NULL)

Arguments

indata

A data.frame object containing the data to be modelled.

dtafile

The file name of the dataset to be imported into MLwiN, which is in Stata format (i.e. with extension .dta).

oldsyntax

Specified as FALSE if new syntax has been used in Formula object, else specified as TRUE (to enable backcompatibility).

resp

A character string (vector) of response variable(s).

levID

A character string (vector) of the specified level ID(s). The ID(s) should be sorted in the descending order of levels (e.g. levID = c('level2', 'level1') where 'level2' is the higher level).

expl

A character string (vector) of explanatory (predictor) variable(s).

rp

A character string (vector) of random part of random variable(s).

D

A character string/vector specifying the type of distribution to be modelled, which can include 'Normal' (the default), 'Binomial', 'Poisson', 'Unordered Multinomial', 'Ordered Multinomial', 'Multivariate Normal', or 'Mixed'. In the case of the latter, 'Mixed' precedes the response types which also need to be listed in D, e.g. c('Mixed', 'Normal', 'Binomial'); these need to be be listed in the same order to which they are referred to in the Formula object (see runMLwiN, Formula.translate, Formula.translate.compat. 'Mixed' combinations can only consist of 'Normal' and 'Binomial' for MCMC estimation.

nonlinear

A character vector specifying linearisation method for IGLS starting values for discrete response models (see Chapter 9 of Rasbash et al 2012, and Goldstein 2011). N = 0 specifies marginal quasi-likelihood linearization (MQL), whilst N = 1 specifies penalised quasi- likelihood linearization (PQL); M = 1 specifies first order approximation, whilst M = 2 specifies second order approximation. nonlinear = c(N = 0, M = 1) by default. First order marginal quasi-likelihood (MQL1) only option for single-level discrete response models.

categ

Specifies categorical variable(s) as a matrix. Each column corresponds to a categorical variable; the first row specifies the name(s) of variable(s); the second row specifies the name(s) of reference group(s), NA(s) if no reference group; the third row states the number of categories for each variable.

notation

Specifies the model subscript notation to be used in the MLwiN equations window. 'class' means no multiple subscripts, whereas 'level' has multiple subscripts.

nonfp

Removes the fixed part of random variable(s). NA if no variable is removed.

clre

A matrix used to estimate some, but not all, of the variances and covariances for a set of coefficients at a particular level. Remove from the random part at level <first row> the covariance matrix element(s) defined by the pair(s) of rows <second row> <third row>. Each row corresponds to a removed entry of the covariance matrix.

Meth

Specifies the maximum likelihood estimation method to be used when generating starting values via (R)IGLS. If Meth = 0 estimation method is set to RIGLS. If Meth = 1 estimation method is set to IGLS (the default setting). If Meth is absent, alternate between IGLS and RIGLS.

merr

A vector which sets-up measurement errors on predictor variables. The first element N defines the number of variables that have measurement errors. Then, for each variable with measurement error, a pair of inputs is required: the first of which is the explanatory variable name as a character string, and the second of which is the variance of the measurement error for this variable. See demo(MCMCGuide14) for an example.

carcentre

If CAR model (i.e. if car is non-NULL), carcentre = TRUE mean-centres all random effects at that level.

maxiter

When generating starting values via (R)IGLS, a numeric value specifying the total number of iterations, from the start, before IGLS estimation halts (if startval = NULL).

convtol

When generating starting values via (R)IGLS, a numeric value specifying the IGLS convergence criterion, as specified in the tol option within estoptions, where startval = NULL) (see runMLwiN). If value of convtol is m, estimation will be deemed to have converged when the relative change in the estimate for any parameter from one iteration to the next is less than 10(-m). Defaults to value of 2 for m if not otherwise specified.

seed

An integer specifying the random seed in MLwiN.

iterations

An integer specifying the number of iterations after burn-in.

burnin

An integer specifying length of the burn-in.

scale

An integer specifying the scale factor of proposed variances; this number will be multiplied by the estimated parameter variance (from IGLS/RIGLS) to give the proposal distribution variance.

thinning

An integer specifying the frequency with which successive values in the Markov chain are stored. By default thinning = 1.

priorParam

A vector specifying the informative priors used, as output from prior2macro.

refresh

An integer specifying how frequently the parameter estimates are refreshed on the screen during iterations; only applies if debugmode = TRUE in estoptions: see runMLwiN.

fixM

Specifies the fixed effect method: 1 for Gibbs Sampling, 2 for univariate MH Sampling and 3 for multivariate MH Sampling.

residM

Specifies the residual method: 1 for Gibbs Sampling, 2 for univariate MH Sampling and 3 for multivariate MH Sampling.

Lev1VarM

Specifies the level 1 variance method: 1 for Gibbs Sampling, 2 for univariate MH Sampling and 3 for multivariate MH Sampling.

OtherVarM

Specifies the variance method for other levels: 1 for Gibbs Sampling, 2 for univariate MH Sampling and 3 for multivariate MH Sampling.

adaption

adaption = 1 indicates adaptation is to be used; 0 otherwise.

priorcode

A vector indicating which default priors are to be used for the variance parameters. It defaults to c(gamma = 1) in which case Gamma priors are used with MLwiN's defaults of Gamma a value (shape) = 0.001 and Gamma b value (scale) = 0.001, although alternative values for shape and scale can be specified in subsequent elements of the vector, e.g. c(gamma = 1, shape = 0.5, scale = 0.2)). Alternatively c(uniform = 1) specifies Uniform priors on the variance scale. To allow for back-compatibility with deprecated syntax used in versions of R2MLwiN prior to 0.8-2, if priorcode is instead specified as an integer, then 1 indicates that Gamma priors are used, whereas 0 indicates that Uniform priors are used. See the section on 'Priors' in the MLwiN help system for more details on the meaning of these priors.

rate

An integer specifying the acceptance rate (as a percentage); this command is ignored if adaption = 0.

tol

An integer specifying tolerance (as a percentage) for the acceptance rate.

lclo

This command toggles on/off the possible forms of complex level 1 variation when using MCMC. lclo = 0 expresses the level 1 variation as a function of the predictors, whereas lclo = 1 expresses the log of the level 1 precision (1/variance) as a function of the predictors.

mcmcOptions

A list of other MCMC options used. See ‘Details’ below.

fact

A list of objects specified for factor analysis. See ‘Details’ below.

xc

Indicates whether model is cross-classified (TRUE) or nested (FALSE). xc = NULL by default (corresponding to FALSE), unless either mm or car are not null, in which case xc = TRUE.

mm

Specifies the structure of a multiple membership model. Can be a list of variable names, a list of vectors, or a matrix (e.g. see df2matrix). In the case of the former, each element of the list corresponds to a level (classification) of the model, in descending order. If a level is not a multiple membership classification, then NA is specified. Otherwise, lists need to be assigned to mmvar and weights, with the former containing columns specifying the classification units, and the latter containing columns specifying the weights. Ignored if EstM = 0, i.e. only applicable to models estimated via MCMC. mm = NULL by default. Supersedes deprecated xclass. E.g. (from demo(MCMCGuide16)) for logearn ~ 1 + age_40 + sex + parttime + (company|1) + (id|1), if company is a multiple membership classification with the variables indicating the classifications in company, company2, company3, company4 and their weights in weight1, weight2, weight3 and weight4 then mm = list(list(mmvar = list('company', 'company2', 'company3', 'company4'), weights = list('weight1', 'weight2', 'weight3', 'weight4')), NA) with the NA, listed last, corresponding to the level 1 identifier (id).

car

A list specifying structure of a conditional autoregressive (CAR) model. Each element of the list corresponds to a level (classification) of the model, in descending order. If a level is not a spatial classification, then NA is specified. Otherwise, lists need to be assigned to carvar and weights, with the former containing columns specifying the spatial classification units, and the latter containing columns specifying the weights. See demo(MCMCGuide17) for examples. car = NULL by default.

BUGO

If non-NULL uses BUGS for MCMC estimation using files specified in modelfile, initfile and datafile.

mem.init

A vector which sets and displays worksheet capacities for the current MLwiN session according to the value(s) specified. By default, the number of levels is nlev+1; worksheet size in thousands of cells is 6000; the number of columns is 2500; the number of explanatory variables is num_vars+10; the number of group labels is 20. nlev is the number of levels specified by levID, and num_vars is approximately the number of explanatory variables calculated initially.

optimat

This option instructs MLwiN to limit the maximum matrix size that can be allocated by the (R)IGLS algorithm. Specify optimat = TRUE if MLwiN gives the following error message 'Overflow allocating smatrix'. This error message arises if one more higher-level units is extremely large (contains more than 800 lower-level units). In this situation runmlwin's default behaviour is to instruct MLwiN to allocate a larger matrix size to the (R)IGLS algorithm than is currently possible. Specifying optimat = TRUE caps the maximum matrix size at 800 lower-level units, circumventing the MLwiN error message, and allowing most MLwiN functionality.

modelfile

A file name where the BUGS model will be saved in .txt format.

initfile

A file name where the BUGS initial values will be saved in .txt format.

datafile

A file name where the BUGS data will be saved in .txt format.

macrofile

A file name where the MLwiN macro file will be saved.

IGLSfile

A file name where the IGLS estimates will be saved.

MCMCfile

A file name where the MCMC estimates will be saved.

chainfile

A file name where the MCMC chains will be saved.

MIfile

A file name where the missing values will be saved.

resifile

A file name where the residual estimates will be saved.

resi.store

A logical value to indicate if residuals are to be stored (TRUE) or not (FALSE).

resioptions

A string vector to specify the various residual options. The 'variance' option calculates the posterior variances instead of the posterior standard errors; the 'standardised' option calculates standardised residuals.

resichains

A file name where the residual chains will be saved.

FACTchainfile

A file name where the factor chains will be saved.

resi.store.levs

An integer vector indicating the levels at which the residual chains are to be stored.

debugmode

A logical value determining whether MLwiN is run in the background or not. The default value is FALSE: i.e. MLwiN is run in the background. If TRUE the MLwiN GUI is opened, and then pauses after the model has been set-up, allowing user to check starting values; pressing 'Resume macro' will then fit the model. Once fit, pressing 'Resume macro' once more will save the outputs to the workdir ready to be read by R2MLwiN. Users can instead opt to 'Abort macro' in which case the outputs are not saved to the workdir. This option currently works for 32 bit version of MLwiN only (automatically switches unless MLwiNPath or options(MLwiNPath) has been set directly to the executable).

startval

A list of numeric vectors specifying the starting values when using MCMC. FP.b corresponds to the estimates for the fixed part; FP.v specifies the variance/covariance estimates for the fixed part; RP.b specifies the variance estimates for the random part; RP.v corresponds to the variance/covariance matrix of the variance estimates for the random part.

dami

This command outputs a complete (i.e. including non-missing responses) response variable y. If dami = c(0, <iter1>, <iter2>,...) then the response variables returned will be the value of y at the iterations quoted (as integers <iter1>, <iter2>, etc.); these can be used for multiple imputation. If dami = 1 the value of y will be the mean estimate from the iterations produced. dami = 2 is as for dami = 1 but with the standard errors of the estimate additionally being stored.

namemap

A mapping of column names to DTA friendly shorter names

saveworksheet

A list of file names (one for each chain) used to store the MLwiN worksheet after the model has been estimated.

Details

A list of other MCMC options as used in the argument mcmcOptions:

  • orth: If orth = 1, orthogonal fixed effect vectors are used; zero otherwise.

  • hcen: An integer specifying the level where we use hierarchical centering.

  • smcm: If smcm = 1, structured MCMC is used; zero otherwise.

  • smvn: If smvn = 1, the structured MVN framework is used; zero otherwise.

  • paex: A matrix of Nx2; in each row, if the second digit is 1, parameter expansion is used at level <the first digit>.

  • mcco: This command allows the user to have constrained settings for the lowest level variance matrix in a multivariate Normal model. If value is 0, it estimates distinct variances for each residual error and distinct covariances for each residual error pair. Four other settings are currently available:

    1 fits stuctured errors with a common correlation paramater and a common variance parameter;
    2 fits AR1 errors with a common variance parameter;
    3 fits structured errors with a common correlation parameter and independent variance parameters;
    4 fits AR1 errors with independent variance parameters.

A list of objects specified for cross-classified and/or multiple membership models, as used in the argument xclass:

  • class: An integer (vector) of the specified class(es).

  • N1: This defines a multiple membership across N1 units at level class. N1>1 if there is multiple membership.

  • weight: If there is multiple membership then the column number weight, which is the length of the dataset, will contain the first set of weights for the multiple membership. Note that there should be N1 weight columns and they should be sequential in the worksheet starting from weight.

  • id: If the response is multivariate then the column number id must be input and this contains the first set of identifiers for the classification. Note that for a p-variate model each lowest level unit contains p records and the identifiers (sequence numbers) for each response variate need to be extracted into id and following columns. There should be N1 of these identifier columns and they should be sequential starting from id in the multivariate case.

  • car: car = TRUE indicates the spatial CAR model; FALSE otherwise. car = FALSE if ignored.

A list of objects specified for factor analysis, as used in the argument fact:

  • nfact: Specifies the number of factors

  • lev.fact: Specifies the level/classification for the random part of the factor for each factor.

  • nfactcor: Specifies the number of correlated factors

  • factcor: a vector specifying the correlated factors: the first element corresponds to the first factor number, the second to the second factor number, the third element corresponds to the starting value for the covariance and the fourth element to whether this covariance is constrained (1) or not (0). If more than one pair of factors is correlated, then repeat this sequence for each pair.

  • loading: A matrix specifying the starting values for the factor loadings and the starting value of the factor variance. Each row corresponds to a factor.

  • constr: A matrix specifying indicators of whether the factor loadings and the factor variance are constrained (1) or not (0).

Value

Outputs a modified version of namemap containing newly generated short names.

Note

Note that for FixM, residM, Lev1VarM and

OtherVarM, not all combinations of methods are available for all sets of parameters and all models.

Author(s)

Zhang, Z., Charlton, C.M.J., Parker, R.M.A., Leckie, G., and Browne, W.J. (2016) Centre for Multilevel Modelling, University of Bristol.

References

Goldstein, H. (2011) Multilevel Statistical Models. 4th Edition. London: John Wiley and Sons.

Rasbash, J., Steele, F., Browne, W.J. and Goldstein, H. (2012) A User's Guide to MLwiN Version 2.26. Centre for Multilevel Modelling, University of Bristol.

See Also

write.IGLS

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.