Nothing
R/write.MCMC.R
In R2MLwiN: Running 'MLwiN' from Within R
Defines functions
write.MCMC
Documented in
write.MCMC
#' Writes MLwiN macros to fit models using Markov chain Monte Carlo (MCMC)
#' methods
#'
#' write.MCMC is an internal function which creates an MLwiN macro file to
#' fit models using MCMC.
#'
#' @param indata A data.frame object containing the data to be modelled.
#' @param dtafile The file name of the dataset to be imported into MLwiN, which
#' is in Stata format (i.e. with extension .dta).
#' @param oldsyntax Specified as \code{FALSE} if new syntax has been used in
#' \code{Formula} object, else specified as \code{TRUE} (to enable
#' backcompatibility).
#' @param resp A character string (vector) of response variable(s).
#' @param 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.
#' \code{levID = c('level2', 'level1')} where \code{'level2'} is the higher
#' level).
#' @param expl A character string (vector) of explanatory (predictor)
#' variable(s).
#' @param rp A character string (vector) of random part of random variable(s).
#' @param D A character string/vector specifying the type of distribution to be modelled, which
#' can include \code{'Normal'} (the default), \code{'Binomial'}, \code{'Poisson'},
#' \code{'Unordered Multinomial'}, \code{'Ordered Multinomial'},
#' \code{'Multivariate Normal'}, or \code{'Mixed'}. In the case of the latter,
#' \code{'Mixed'} precedes the response types which also need to be listed in
#' \code{D}, e.g. \code{c('Mixed', 'Normal', 'Binomial')}; these need to be
#' be listed in the same order to which they are referred to in the
#' \code{Formula} object (see \code{\link{runMLwiN}}, \code{\link{Formula.translate}},
#' \code{\link{Formula.translate.compat}}. \code{'Mixed'} combinations can only consist of
#' \code{'Normal'} and \code{'Binomial'} for MCMC estimation.
#' @param 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).
#' \code{N = 0} specifies marginal quasi-likelihood
#' linearization (MQL), whilst \code{N = 1} specifies penalised quasi-
#' likelihood linearization (PQL); \code{M = 1} specifies first order
#' approximation, whilst \code{M = 2} specifies second order approximation.
#' \code{nonlinear = c(N = 0, M = 1)} by default. First order marginal
#' quasi-likelihood (MQL1) only option for single-level discrete response
#' models.
#' @param 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),
#' \code{NA}(s) if no reference group; the third row states the number of
#' categories for each variable.
#' @param notation Specifies the model subscript notation to be used in the
#' MLwiN equations window. \code{'class'} means no multiple subscripts, whereas
#' \code{'level'} has multiple subscripts.
#' @param nonfp Removes the fixed part of random variable(s). \code{NA} if no
#' variable is removed.
#' @param 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.
#' @param Meth Specifies the maximum likelihood estimation method to be used
#' when generating starting values via (R)IGLS.
#' If \code{Meth = 0} estimation method is set to RIGLS. If \code{Meth = 1}
#' estimation method is set to IGLS (the default setting). If \code{Meth} is
#' absent, alternate between IGLS and RIGLS.
#' @param merr A vector which sets-up measurement errors on predictor
#' variables. The first element \code{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 \code{demo(MCMCGuide14)} for an
#' example.
#' @param carcentre If CAR model (i.e. if \code{car} is non-\code{NULL}),
#' \code{carcentre = TRUE} mean-centres all random effects at that level.
#' @param 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 \code{startval = NULL}).
#' @param convtol When generating starting values via (R)IGLS, a numeric
#' value specifying the IGLS convergence criterion, as
#' specified in the \code{tol} option within \code{estoptions}, where
#' \code{startval = NULL}) (see \code{\link{runMLwiN}}). If value of
#' \code{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 \code{2} for m if not
#' otherwise specified.
#' @param seed An integer specifying the random seed in MLwiN.
#' @param iterations An integer specifying the number of iterations after
#' burn-in.
#' @param burnin An integer specifying length of the burn-in.
#' @param 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.
#' @param thinning An integer specifying the frequency with which successive
#' values in the Markov chain are stored. By default \code{thinning = 1}.
#' @param priorParam A vector specifying the informative priors used, as output
#' from \code{\link{prior2macro}}.
#' @param refresh An integer specifying how frequently the parameter estimates
#' are refreshed on the screen during iterations; only applies if
#' \code{debugmode = TRUE} in \code{estoptions}:
#' see \code{\link[R2MLwiN]{runMLwiN}}.
#' @param fixM Specifies the fixed effect method: \code{1} for Gibbs Sampling,
#' \code{2} for univariate MH Sampling and \code{3} for multivariate MH
#' Sampling.
#' @param residM Specifies the residual method: \code{1} for Gibbs Sampling,
#' \code{2} for univariate MH Sampling and \code{3} for multivariate MH
#' Sampling.
#' @param Lev1VarM Specifies the level 1 variance method: \code{1} for Gibbs
#' Sampling, \code{2} for univariate MH Sampling and \code{3} for multivariate
#' MH Sampling.
#' @param OtherVarM Specifies the variance method for other levels: \code{1}
#' for Gibbs Sampling, \code{2} for univariate MH Sampling and \code{3} for
#' multivariate MH Sampling.
#' @param adaption \code{adaption = 1} indicates adaptation is to be used;
#' \code{0} otherwise.
#' @param priorcode A vector indicating which default priors are to be used
#' for the variance parameters. It defaults to \code{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. \code{c(gamma = 1, shape = 0.5, scale = 0.2)}). Alternatively
#' \code{c(uniform = 1)} specifies Uniform priors on the variance scale. To allow
#' for back-compatibility with deprecated syntax used in versions of
#' \pkg{R2MLwiN} prior to 0.8-2, if \code{priorcode} is instead specified as
#' an integer, then \code{1} indicates that Gamma priors are used, whereas
#' \code{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.
#' @param rate An integer specifying the acceptance rate (as a percentage);
#' this command is ignored if \code{adaption = 0}.
#' @param tol An integer specifying tolerance (as a percentage) for the acceptance rate.
#' @param lclo This command toggles on/off the possible forms of complex level
#' 1 variation when using MCMC. \code{lclo = 0} expresses the level
#' 1 variation as a function of the predictors, whereas \code{lclo = 1} expresses the
#' log of the level 1 precision (1/variance) as a function of the predictors.
#' @param mcmcOptions A list of other MCMC options used. See `Details' below.
#' @param fact A list of objects specified for factor analysis. See `Details'
#' below.
#' @param xc Indicates whether model is cross-classified (\code{TRUE}) or
#' nested (\code{FALSE}). \code{xc = NULL} by default (corresponding to
#' \code{FALSE}), unless either \code{mm} or \code{car} are not null, in
#' which case \code{xc = TRUE}.
#' @param 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
#' \code{\link{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 \code{NA} is specified. Otherwise, lists need to be assigned to
#' \code{mmvar} and \code{weights}, with the former containing columns
#' specifying the classification units, and the latter containing columns
#' specifying the weights. Ignored if \code{EstM = 0}, i.e. only applicable to models estimated via
#' MCMC. \code{mm = NULL} by default. Supersedes deprecated \code{xclass}.
#' E.g. (from \code{demo(MCMCGuide16)}) for
#' \code{logearn ~ 1 + age_40 + sex + parttime + (company|1) + (id|1)}, if
#' \code{company} is a multiple membership classification with the variables
#' indicating the classifications in \code{company}, \code{company2},
#' \code{company3}, \code{company4} and their weights in \code{weight1}, \code{weight2},
#' \code{weight3} and \code{weight4} then
#' \code{mm = list(list(mmvar = list('company', 'company2', 'company3', 'company4'),}
#' \code{weights = list('weight1', 'weight2', 'weight3', 'weight4')), NA)}
#' with the \code{NA}, listed last, corresponding to the level 1 identifier (\code{id}).
#' @param 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 \code{NA} is specified. Otherwise, lists need to be assigned to
#' \code{carvar} and \code{weights}, with the former containing columns
#' specifying the spatial classification units, and the latter containing
#' columns specifying the weights. See \code{demo(MCMCGuide17)} for examples.
#' \code{car = NULL} by default.
#' @param BUGO If non-\code{NULL} uses BUGS for MCMC estimation using files
#' specified in \code{modelfile}, \code{initfile} and \code{datafile}.
#' @param 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 \code{nlev}+1; worksheet size in thousands of cells
#' is 6000; the number of columns is 2500; the number of explanatory variables
#' is \code{num_vars}+10; the number of group labels is 20. \code{nlev} is the
#' number of levels specified by \code{levID}, and \code{num_vars} is
#' approximately the number of explanatory variables calculated initially.
#' @param optimat This option instructs MLwiN to limit the maximum matrix size
#' that can be allocated by the (R)IGLS algorithm. Specify \code{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
#' \code{optimat = TRUE} caps the maximum matrix size at 800 lower-level units,
#' circumventing the MLwiN error message, and allowing most MLwiN
#' functionality.
#' @param modelfile A file name where the BUGS model will be saved in .txt
#' format.
#' @param initfile A file name where the BUGS initial values will be saved
#' in .txt format.
#' @param datafile A file name where the BUGS data will be saved in .txt
#' format.
#' @param macrofile A file name where the MLwiN macro file will be saved.
#' @param IGLSfile A file name where the IGLS estimates will be saved.
#' @param MCMCfile A file name where the MCMC estimates will be saved.
#' @param chainfile A file name where the MCMC chains will be saved.
#' @param MIfile A file name where the missing values will be saved.
#' @param resifile A file name where the residual estimates will be saved.
#' @param resi.store A logical value to indicate if residuals are to be stored
#' (\code{TRUE}) or not (\code{FALSE}).
#' @param resioptions A string vector to specify the various residual options.
#' The \code{'variance'} option calculates the posterior variances instead of
#' the posterior standard errors; the \code{'standardised'} option calculates standardised
#' residuals.
#' @param resichains A file name where the residual chains will be saved.
#' @param FACTchainfile A file name where the factor chains will be saved.
#' @param resi.store.levs An integer vector indicating the levels at which the
#' residual chains are to be stored.
#' @param debugmode A logical value determining whether MLwiN is run in the
#' background or not. The default value is \code{FALSE}: i.e. MLwiN is run in
#' the background. If \code{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 \code{workdir} ready to be read by \pkg{R2MLwiN}. Users can
#' instead opt to 'Abort macro' in which case the outputs are not saved to the
#' \code{workdir}. This option currently
#' works for 32 bit version of MLwiN only (automatically switches unless
#' \code{MLwiNPath} or \code{options(MLwiNPath)}
#' has been set directly to the executable).
#' @param startval A list of numeric vectors specifying the starting values
#' when using MCMC. \code{FP.b} corresponds to the estimates for the fixed
#' part; \code{FP.v} specifies the variance/covariance estimates for the fixed
#' part; \code{RP.b} specifies the variance estimates for the random part;
#' \code{RP.v} corresponds to the variance/covariance matrix of the variance
#' estimates for the random part.
#' @param dami This command outputs a complete (i.e. including non-missing
#' responses) response variable y. If \code{dami = c(0, <iter1>, <iter2>,...)} then
#' the response variables returned will be the value of y at the iterations
#' quoted (as integers \code{<iter1>, <iter2>}, etc.); these can be used for
#' multiple imputation. If \code{dami = 1} the value of y will be the mean
#' estimate from the iterations produced. \code{dami = 2} is as for \code{dami = 1}
#' but with the standard errors of the estimate additionally being stored.
#' @param namemap A mapping of column names to DTA friendly shorter names
#' @param saveworksheet A list of file names (one for each chain) used to store the
#' MLwiN worksheet after the model has been estimated.
#' @param rng.version An integer indicating which random number generator should
#' be used by MLwiN.
#'
#' @details A list of other MCMC options as used in the argument
#' \code{mcmcOptions}:
#' \itemize{
#' \item \code{orth}: If \code{orth = 1}, orthogonal fixed effect
#' vectors are used; zero otherwise.
#' \item \code{hcen}: An integer specifying the
#' level where we use hierarchical centering.
#' \item \code{smcm}: If \code{smcm = 1},
#' structured MCMC is used; zero otherwise.
#' \item \code{smvn}: If \code{smvn = 1}, the
#' structured MVN framework is used; zero otherwise.
#' \item \code{paex}: A matrix of Nx2; in each row, if the second digit is \code{1}, parameter expansion
#' is used at level <the first digit>.
#' \item \code{mcco}: This
#' command allows the user to have constrained settings for the lowest level
#' variance matrix in a multivariate Normal model. If value is \code{0},
#' it estimates distinct variances for each residual error and distinct covariances
#' for each residual error pair. Four other
#' settings are currently available:\cr
#' \tabular{ll}{\code{1} \tab fits stuctured errors with a common correlation paramater and a common variance parameter;\cr
#' \code{2} \tab fits AR1 errors with a common variance parameter;\cr \code{3} \tab fits structured errors with a common
#' correlation parameter and independent variance parameters;\cr \code{4} \tab fits AR1 errors with independent variance
#' parameters.\cr }
#' }
#'
#' A list of objects specified for cross-classified and/or multiple membership
#' models, as used in the argument \code{xclass}:
#' \itemize{
#' \item \code{class}: An integer
#' (vector) of the specified class(es).
#' \item \code{N1}: This defines a multiple
#' membership across \code{N1} units at level \code{class}. \code{N1}>1 if
#' there is multiple membership.
#' \item \code{weight}: If there is multiple
#' membership then the column number \code{weight}, which is the length of the
#' dataset, will contain the first set of weights for the multiple membership.
#' Note that there should be \code{N1} weight columns and they should be
#' sequential in the worksheet starting from \code{weight}.
#' \item \code{id}: If the
#' response is multivariate then the column number \code{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 \code{id} and following columns. There should be \code{N1} of
#' these identifier columns and they should be sequential starting from
#' \code{id} in the multivariate case.
#' \item \code{car}: \code{car = TRUE} indicates
#' the spatial CAR model; \code{FALSE} otherwise. \code{car = FALSE} if ignored.
#' }
#'
#' A list of objects specified for factor analysis, as used in the argument
#' \code{fact}:
#' \itemize{
#' \item \code{nfact}: Specifies the number of factors
#' \item \code{lev.fact}: Specifies the level/classification for the random part of
#' the factor for each factor.
#' \item \code{nfactcor}: Specifies the number of
#' correlated factors
#' \item \code{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
#' (\code{1}) or not (\code{0}). If more than one pair of factors is correlated,
#' then repeat this sequence for each pair.
#' \item \code{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.
#' \item \code{constr}: A matrix
#' specifying indicators of whether the factor loadings and the factor variance
#' are constrained (\code{1}) or not (\code{0}).
#' }
#'
#' @return Outputs a modified version of namemap containing newly generated
#' short names.
#' @note Note that for \code{FixM}, \code{residM}, \code{Lev1VarM} and
#'
#' \code{OtherVarM}, not all combinations of methods are available for all sets
#' of parameters and all models.
#'
#' @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.
#'
#' @author Zhang, Z., Charlton, C.M.J., Parker, R.M.A., Leckie, G., and Browne,
#' W.J. (2016) Centre for Multilevel Modelling, University of Bristol.
#'
#' @seealso \code{\link{write.IGLS}}
#'
write.MCMC <- function(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, rng.version = 10) {
shortname <- function(...) {
name <- paste0(...)
if (!name %in% names(namemap)) {
sname <- paste0("v", digest::digest(name, algo = "xxhash64", serialize = FALSE))
names(sname) <- name
namemap <<- c(namemap, sname)
}
return(namemap[[name]])
}
nlev <- length(levID)
nrp <- length(rp)
if (nrp > 0) {
rp.names <- names(rp)
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
for (i in 1:nrp) {
temp <- rp.names[i]
rp.names[i] <- paste("rp", as.integer(sub("rp", "", temp)) + 1, sep = "")
}
}
names(rp) <- rp.names
}
num.expl.init <- function(p, nonfp, categ) {
num_vars <- 0
if (is.na(nonfp[1])) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
num_vars <- num_vars + 1
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])])
} else {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1
}
}
} else {
if (sum(p == nonfp) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
num_vars <- num_vars + 1
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])])
} else {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1
}
}
}
}
num_vars
}
if (is.list(expl)) {
sep.coeff <- expl$sep.coeff
if (is.list(nonfp)) {
nonfp.sep <- nonfp$nonfp.sep
nonfp.common <- nonfp$nonfp.common
}
if (!is.na(sep.coeff[1])) {
num_vars <- sum(sapply(sep.coeff, function(x) num.expl.init(x, nonfp.sep, categ)))
} else {
num_vars <- 0
}
if (D[1] == "Multinomial") {
nresp <- length(levels(indata[, resp])) - 1
## it is true if adding each expl variables separately
num_vars <- num_vars * nresp
}
if (D[1] == "Multivariate Normal") {
nresp <- length(resp)
num_vars <- num_vars * nresp
}
common.coeff <- expl$common.coeff
num_vars <- num_vars + sum(sapply(common.coeff, function(x) num.expl.init(x, nonfp$nonfp.common, categ)))
common.coeff.id <- expl$common.coeff.id
} else {
num_vars <- sum(sapply(expl, function(x) num.expl.init(x, nonfp, categ)))
if (D[1] == "Multinomial") {
nresp <- length(levels(indata[, resp])) - 1
num_vars <- num_vars * nresp
}
if (D[1] == "Multivariate Normal") {
nresp <- length(resp)
num_vars <- num_vars * nresp
}
}
if (nrp > 0) {
for (ii in 1:nrp) {
rpx <- rp[[rp.names[ii]]]
nrpx <- length(rpx)
if (nrpx == 1)
num_vars <- num_vars + 1
if (nrpx >= 2)
num_vars <- num_vars + nrpx * (nrpx - 1)/2 + nrpx
}
}
## Write into macro file
wrt <- function(x) write(x, macrofile, append = TRUE)
cat(file = macrofile)
wrt("ECHO 0")
wrt("NOTE *****************************************************************")
wrt("NOTE MLwiN macro created by rtomlwin command")
wrt("NOTE *****************************************************************")
wrt("")
wrt("NOTE Initialise MLwiN storage")
if (mem.init[1] == "default") {
wrt(paste("INIT ", nlev + 1, " 6000 2500 ", num_vars + 10, " 30", sep = ""))
} else {
wrt(paste("INIT ", mem.init[1], mem.init[2], mem.init[3], mem.init[4], mem.init[5]))
}
wrt("NOTE Limit the maximum matrix size")
if (optimat) {
wrt("OPTS 1")
} else {
wrt("OPTS 0")
}
wrt("MARK 0")
wrt("NOTE Import the R data into MLwiN")
wrt(paste("RSTA '", dtafile, "'", sep = ""))
wrt("NOTE Set the random number generator version")
wrt(paste0("RNGV ", rng.version))
wrt("NOTE Correct column names")
for (name in names(namemap)) {
wrt(paste0("COLN '", namemap[[name]], "' b50"))
wrt(paste0("DESC cb50 '", name, "'"))
wrt(paste0("NAME cb50 '", name, "'"))
}
if ("class" %in% notation) {
wrt("INDE 1")
}
if (!(D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed")) {
wrt("NOTE Specify the response variable")
for (p in resp) wrt(paste("RESP '", p, "'", sep = ""))
wrt("")
}
if (D[1] == "Binomial") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("RDISt 1 0")
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
if (D[2] == "logit") {
wrt("LFUN 0")
DD2 <- 0
}
if (D[2] == "probit") {
wrt("LFUN 1")
DD2 <- 1
}
if (D[2] == "cloglog") {
wrt("LFUN 2")
DD2 <- 2
}
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (D[1] == "Mixed") {
nresp <- length(resp)
for (ii in 1:nresp) wrt(paste("MVAR 1 '", resp[ii], "'", sep = ""))
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:2
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("IDEN 1 'resp_indicator'")
jj <- 1
for (ii in 2:length(D)) {
if (D[[ii]][1] == "Binomial") {
wrt(paste("RDISt ", jj, " 0", sep = ""))
wrt(paste("DOFFs ", jj, " '", D[[ii]][3], "'", sep = ""))
if (D[[ii]][2] == "logit") {
wrt("LFUN 0")
DD2 <- 0
}
if (D[[ii]][2] == "probit") {
wrt("LFUN 1")
DD2 <- 1
}
if (D[[ii]][2] == "cloglog") {
wrt("LFUN 2")
DD2 <- 2
}
}
if (D[[ii]][1] == "Poisson") {
wrt(paste("RDISt ", jj, " 1", sep = ""))
wrt("LFUN 3")
DD2 <- 3
DD2 <- 3
if (!is.na(D[[jj]][3])) {
wrt(paste("DOFFs ", jj, " '", D[[jj]][3], "'", sep = ""))
}
}
jj <- jj + 1
}
wrt("")
if (is.list(expl)) {
if (!is.na(sep.coeff[1])) {
interpos1 <- grep("\\:", sep.coeff)
if (!isTRUE(oldsyntax) || length(interpos1) == 0) {
for (x in 1:length(sep.coeff)) {
p <- sep.coeff[x]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- sep.coeff[interpos1]
explx <- sep.coeff[-interpos1]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
sep.coeff <- c(explx, exply)
}
}
interpos2 <- grep("\\:", common.coeff)
if (!isTRUE(oldsyntax) || length(interpos2) == 0) {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (!(y %in% interpos2)) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
} else {
TT <- ""
interx <- unlist(strsplit(p, "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
wrt("RPAT")
common.coeff[y] <- paste(p, ".", bb, sep = "")
}
}
} else {
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
}
wrt("")
}
if (D[1] == "Multivariate Normal") {
nresp <- length(resp)
for (ii in 1:nresp) wrt(paste("MVAR 1 '", resp[ii], "'", sep = ""))
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:2
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("IDEN 1 'resp_indicator'")
wrt("LFUN 0")
wrt("")
if (is.list(expl)) {
if (!is.na(sep.coeff[1])) {
interpos1 <- grep("\\:", sep.coeff)
if (!isTRUE(oldsyntax) || length(interpos1) == 0) {
for (x in 1:length(sep.coeff)) {
p <- sep.coeff[x]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- sep.coeff[interpos1]
explx <- sep.coeff[-interpos1]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
sep.coeff <- c(explx, exply)
}
}
interpos2 <- grep("\\:", common.coeff)
if (!isTRUE(oldsyntax) || length(interpos2) == 0) {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (!(y %in% interpos2)) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
} else {
TT <- ""
interx <- unlist(strsplit(p, "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
wrt("RPAT")
common.coeff[y] <- paste(p, ".", bb, sep = "")
}
}
} else {
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
}
wrt("")
}
if (D[1] == "Multinomial") {
wrt("LINK 2 G21")
wrt("NAME G21[1] 'resp' G21[2] 'resp_indicator'")
wrt("LINK 0 G21")
wrt(paste("MNOM ", as.numeric(D[4]), " '", resp, "' ", "'resp' 'resp_indicator' ", as.numeric(D[5]), sep = ""))
wrt("RESP 'resp'")
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:c(nlev - 1)) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("IDEN 1 'resp_indicator'")
wrt("")
if (as.numeric(D[4]) == 0)
wrt("RDISt 1 4") else wrt("RDISt 1 5")
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
if (D[2] == "logit") {
wrt("LFUN 0")
DD2 <- 0
}
if (D[2] == "probit") {
wrt("LFUN 1")
DD2 <- 1
}
if (D[2] == "cloglog") {
wrt("LFUN 2")
DD2 <- 2
}
wrt("")
if (is.list(expl)) {
if (!is.na(sep.coeff[1])) {
interpos1 <- grep("\\:", sep.coeff)
if (!isTRUE(oldsyntax) || length(interpos1) == 0) {
for (x in 1:length(sep.coeff)) {
p <- sep.coeff[x]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- sep.coeff[interpos1]
explx <- sep.coeff[-interpos1]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
sep.coeff <- c(explx, exply)
}
}
interpos2 <- grep("\\:", common.coeff)
if (!isTRUE(oldsyntax) || length(interpos2) == 0) {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (!(y %in% interpos2)) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
} else {
TT <- ""
interx <- unlist(strsplit(p, "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
wrt("RPAT")
common.coeff[y] <- paste(p, ".", bb, sep = "")
}
}
} else {
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
}
wrt("")
}
if (D[1] == "Normal") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("NOTE Specify covariate(s) used anywhere in the model")
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (D[1] == "Poisson") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("RDISt 1 1")
wrt("LFUN 3")
DD2 <- 3
if (!is.na(D[3])) {
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
}
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (D[1] == "Negbinom") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("RDISt 1 2")
wrt("LFUN 3")
DD2 <- 3
if (!is.na(D[3])) {
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
}
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (is.list(nonfp)) {
wrt("NOTE Turn off the fixed part of the explanatary variable(s)")
nonfp.sep <- nonfp$nonfp.sep
nonfp.common <- nonfp$nonfp.common
if (!is.na(nonfp.sep[1])) {
interpos <- grep("\\:", nonfp.sep)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in nonfp.sep) wrt(paste("FPAR 0 '", p, "'", sep = ""))
} else {
for (i in 1:length(nonfp.sep)) {
if (i %in% interpos && oldsyntax) {
wrt(paste("FPAR 0 '", gsub("\\:", "\\.", nonfp.sep[i]), "'", sep = ""))
} else {
wrt(paste("FPAR 0 '", nonfp.sep[i], "'", sep = ""))
}
}
}
}
if (!is.na(nonfp.common[1])) {
interpos <- grep("\\:", nonfp.common)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in nonfp.common) wrt(paste("FPAR 0 '", p, "'", sep = ""))
} else {
for (i in 1:length(nonfp.common)) {
if (i %in% interpos && oldsyntax) {
wrt(paste("FPAR 0 '", gsub("\\:", "\\.", nonfp.common[i]), "'", sep = ""))
} else {
wrt(paste("FPAR 0 '", nonfp.common[i], "'", sep = ""))
}
}
}
}
} else {
if (!is.na(nonfp[1])) {
wrt("NOTE Turn off the fixed part of the explanatory varible(s)")
for (p in nonfp) {
if(oldsyntax){
p <- gsub("\\:", "\\.", p)
}
wrt(paste("FPAR 0 '", p, "'", sep = ""))
}
}
}
wrt("")
wrt("NOTE Specify random part covariate(s)")
if (nrp > 0) {
for (ii in 1:nrp) {
for (p in rp[[ii]]) {
if (oldsyntax) {
p <- gsub("\\:", "\\.", p)
}
wrt(paste("SETV ", as.numeric(sub("rp", "", rp.names[ii])), " '", p, "'", sep = ""))
}
}
}
if (!is.null(clre)) {
nclre <- ncol(clre)
for (ii in 1:nclre) {
if (oldsyntax){
wrt(paste("CLRE ", as.numeric(clre[1, ii]), " '", gsub("\\:", "\\.", clre[2, ii]), "' '",
gsub("\\:", "\\.", clre[3, ii]), "'", sep = ""))
} else {
wrt(paste("CLRE ", as.numeric(clre[1, ii]), " '", clre[2, ii], "' '", clre[3, ii], "'", sep = ""))
}
}
}
nexpl <- length(expl)
wrt("")
wrt("NOTE Set estimation method")
if (Meth != 2) {
wrt(paste("METH", Meth))
}
wrt(paste("LINE ", nonlinear[1], nonlinear[2]))
wrt("")
if (!is.null(fact)) {
TT <- NULL
for (i in 1:fact$nfact) {
TT <- c(TT, fact$lev.fact[i] + 1, matrix(rbind(fact$loading[i, ], fact$constr[i, ]), nrow = 1))
}
if (fact$nfactcor > 0)
TT <- c(TT, fact$factcor)
FACT <- as.vector(c(length(resp), fact$nfact, fact$nfactcor, TT))
rm(TT)
TT <- ""
for (i in 1:length(FACT)) {
TT <- paste(TT, FACT[i])
}
wrt(paste("FACT ", TT))
wrt("LINK 2 G21")
wrt("SMFA 1 G21[1]")
wrt("SMFA 2 G21[2]")
}
if (D[1] == "Normal") {
wrt("PREF 0")
wrt("POST 0")
}
if ("simple" %in% notation) {
wrt("EXISt 'cons' b1000")
wrt("SWITch b1000")
wrt("CASE 0:")
wrt("EXISt 'Intercept' b1000")
wrt("SWITch b1000")
wrt("CASE 1:")
wrt("COLN 'Intercept' b1000")
wrt("NAME cb1000 'cons'")
wrt("ENDSWITch")
wrt("ENDSWITch")
wrt("NOTA 1")
}
wrt("NAME c1098 '_FP_b'")
wrt("NAME c1099 '_FP_v'")
wrt("NAME c1096 '_RP_b'")
wrt("NAME c1097 '_RP_v'")
wrt("NOTE Fit the model")
wrt("ECHO 1")
wrt("BATC 1")
wrt("MAXI 2")
wrt("STAR")
if (!is.null(startval)) {
if (!is.null(startval$FP.b)) {
wrt(paste("JOIN ", paste(startval$FP.b, collapse = " "), " '_FP_b'", sep = ""))
}
if (!is.null(startval$FP.v)) {
wrt(paste("JOIN ", paste(startval$FP.v[!upper.tri(startval$FP.v)], collapse = " "), " '_FP_v'", sep = ""))
}
if (!is.null(startval$RP.b)) {
wrt(paste("JOIN ", paste(startval$RP.b, collapse = " "), " '_RP_b'", sep = ""))
if (D[1] == "Multinomial" && D[4] == 0) {
wrt("JOIN '_RP_b' 1 '_RP_b'")
}
}
if (!is.null(startval$RP.v)) {
wrt(paste("JOIN ", paste(startval$RP.v[!upper.tri(startval$RP.v)], collapse = " "), " '_RP_v'", sep = ""))
if (D[1] == "Multinomial" && D[4] == 0) {
wrt(paste("JOIN '_RP_v' ", paste(rep(0, length(startval$RP.b)+1), collapse = " "), " '_RP_v'"))
}
}
} else {
wrt(paste("TOLE", convtol))
wrt(paste("MAXI", maxiter))
wrt("BATC 1")
wrt("NEXT")
}
wrt("ECHO 0")
wrt("ITNU 0 b21")
wrt("CONV b22")
wrt("")
wrt("NOTE *****************************************************************")
wrt("")
wrt("NOTE *****************************************************************")
wrt("NOTE Export the model results to R")
wrt("NOTE *****************************************************************")
wrt("LINK 1 G30")
wrt("NAME G30[1] '_Stats'")
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
wrt("NOBS 2 b31 b32")
} else {
wrt("NOBS 1 b31 b32")
}
wrt("EDIT 1 '_Stats' b31")
wrt("EDIT 2 '_Stats' b32")
if (D[1] == "Normal" || D[1] == "Multivariate Normal") {
wrt("LIKE b100")
}
wrt("EDIT 3 '_Stats' b100")
wrt("EDIT 7 '_Stats' b21")
wrt("EDIT 8 '_Stats' b22")
wrt("NAME c1094 '_esample'")
wrt("SUM '_esample' b1")
wrt("EDIT 9 '_Stats' b1")
wrt(paste("PSTA '", IGLSfile, "' ", "'_FP_b' ", "'_FP_v' ", "'_RP_b' ", "'_RP_v' ", "'_Stats'", sep = ""))
wrt("LINK 0 G30")
wrt("NOTE *****************************************************************")
wrt("NOTE Set estimation method to MCMC")
wrt("NOTE *****************************************************************")
wrt("EMODe 3")
if (!is.null(resi.store.levs)) {
wrt(paste0("LINK ", length(resi.store.levs), " G22"))
for (i in 1:length(resi.store.levs)) {
resilev <- resi.store.levs[i]
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
resilev <- resilev + 1
}
wrt(paste("SMRE ", resilev, " G22[", i, "]", sep = ""))
wrt(paste("NAME ", " G22[", i, "] 'resi_lev", resi.store.levs[i], "'", sep = ""))
}
}
if (!is.null(merr)) {
nmerr <- as.numeric(merr[1])
tt <- paste("MERR ", nmerr)
j <- 1
for (ii in 1:nmerr) {
tt <- paste(tt, " '", merr[j + 1], "' ", as.numeric(merr[j + 2]), sep = "")
j <- j + 2
}
wrt(tt)
}
wrt(paste("ORTH", mcmcOptions$orth))
if (mcmcOptions$hcen == 0) {
wrt(paste("HCEN", mcmcOptions$hcen))
} else {
wrt(paste("HCEN", 1, mcmcOptions$hcen))
}
wrt(paste("SMCM", mcmcOptions$smcm))
wrt(paste("SMVN", mcmcOptions$smvn))
if (is.matrix(mcmcOptions$paex)) {
apply(mcmcOptions$paex, 1, function(x) wrt(paste("PAEX", x[1], x[2])))
} else {
wrt(paste("PAEX", mcmcOptions$paex[1], mcmcOptions$paex[2]))
}
if (D[1] == "Multivariate Normal")
wrt(paste("MCCO ", mcmcOptions$mcco))
if (!is.null(xc) && isTRUE(xc)) {
wrt("XCLA 1")
}
if (!is.null(mm)) {
for (i in 1:length(mm)) {
if (!any(is.na(mm[[i]]))) {
if (D[1] != "Multivariate Normal") {
mmlev <- (length(mm) - i) + 1
wrt(paste0("MULM ", mmlev, " ", length(mm[[i]]$mmvar), " '", mm[[i]]$weights[[1]], "'"))
} else {
mmlev <- (length(mm) - i) + 2
wrt(paste0("MULM ", mmlev, " ", length(mm[[i]]$mmvar), " '", mm[[i]]$weights[[1]], "' '", mm[[i]]$mmvar[[1]],
"'"))
}
}
}
}
if (!is.null(car)) {
for (i in 1:length(car)) {
if (!any(is.na(car[[i]]))) {
carlev <- (length(car) - i) + 1
wrt(paste0("MULM ", carlev, " ", length(car[[i]]$carvar), " '", car[[i]]$weights[[1]], "' '", car[[i]]$carvar[[1]],
"'"))
wrt(paste("CARP", carlev, "1"))
}
}
if (!is.null(carcentre) && isTRUE(carcentre)) {
wrt("CARC 1")
}
}
wrt("")
wrt("NOTE Set MCMC seed")
wrt(paste("MCRS ", seed, sep = ""))
wrt("")
wrt("NOTE Set prior distribution parameters")
if (priorParam[1] != "default") {
wrt("PRIOR c1092")
lenpp <- length(priorParam)
tempt <- " "
for (i in 1:lenpp) tempt <- paste(tempt, priorParam[i])
wrt(paste("JOIN", tempt, " c1092", sep = ""))
wrt("")
} else {
wrt("")
}
if (D[1] == "Normal") {
wrt("PREF 0")
wrt("POST 0")
}
if (nlev > 1 && !isTRUE(xc)) {
wrt("MISR 0")
wrt("LINK 2 G30")
if (D[1] == "Multinomial" && as.numeric(D[4]) == 0) {
len.rpx <- 2
wrt(paste0("LINK ", len.rpx, " G27"))
wrt("LINK 1 G28")
wrt("NOTE Calculate MCMC starting values for level 2 residuals")
wrt("RLEV 2")
wrt("RFUN")
wrt("RCOV 2")
wrt("ROUT G27 G28")
wrt("RESI")
wrt("JOIN G30[1] G27 G30[1]")
wrt("JOIN G30[2] G28 G30[2]")
wrt("ERAS G27")
wrt("ERAS G28")
wrt("LINK 0 G27")
wrt("LINK 0 G28")
}
if (nrp > 0) {
for (j in nrp:1) {
if (as.numeric(sub("rp", "", rp.names[j])) != 1) {
rpx <- rp[[j]]
len.rpx <- length(rpx)
wrt(paste0("LINK ", len.rpx, " G27"))
wrt("LINK 1 G28")
wrt(paste("NOTE Calculate MCMC starting values for level ", as.numeric(sub("rp", "", rp.names[j])),
" residuals", sep = ""))
wrt(paste("RLEV ", as.numeric(sub("rp", "", rp.names[j])), sep = ""))
wrt("RFUN")
wrt("RCOV 2")
wrt("ROUT G27 G28")
wrt("RESI")
wrt("JOIN G30[1] G27 G30[1]")
wrt("JOIN G30[2] G28 G30[2]")
wrt("ERAS G27")
wrt("ERAS G28")
wrt("LINK 0 G27")
wrt("LINK 0 G28")
}
}
}
wrt("MISR 1")
}
if (D[1] == "Normal")
DD <- 1
if (D[1] == "Binomial")
DD <- 2
if (D[1] == "Mixed")
DD <- 5
if (D[1] == "Poisson")
DD <- 3
if (D[1] == "Multivariate Normal")
DD <- 4
if (D[1] == "Multinomial") {
if (as.numeric(D[4]) == 0)
DD <- 6 else DD <- 7
wrt("CLRV 2")
}
if (D[1] == "Negbinom")
DD <- 8
wrt(paste("LCLO ", lclo, sep = ""))
if (priorcode["gamma"] == 0) {
wrt("UNIP")
} else {
if (length(priorcode) > 1) {
wrt(paste("GAMP", priorcode["shape"], priorcode["scale"]))
}
}
if (debugmode) {
wrt("NOTE Open the equations window")
wrt("WSET 15 1")
wrt("EXPA 3")
wrt("ESTM 2")
wrt("PAUS")
}
priorcol <- ""
if (priorParam[1] != "default") {
priorcol <- "c1092"
}
if ((!is.null(BUGO)) && !(D[1] == "Mixed") && nrp > 0) {
if (D[1] == "Normal" || D[1] == "Multivariate Normal")
DD2 <- 0
if (nlev > 1 && !isTRUE(xc)) {
wrt(paste("BUGO 6 ", DD, " ", DD2, " G30[1] ", priorcol, " '", modelfile, "' ", "'", initfile, "' ", "'",
datafile, "'", sep = ""))
wrt("ERAS G30")
wrt("LINK 0 G30")
} else {
wrt(paste("BUGO 6 ", DD, " ", DD2, " ", "'", modelfile, "' ", "'", initfile, "' ", "'", datafile, "'",
sep = ""))
}
} else {
wrt("NOTE fit the model in MCMC")
wrt(paste("MTOT ", iterations, sep = ""))
wrt("ECHO 1")
residcols <- ""
if (nlev > 1 && !isTRUE(xc)) {
residcols <- "G30[1] G30[2]"
}
wrt(paste("MCMC 0", burnin, adaption, scale, rate, tol, residcols, priorcol, fixM, residM, Lev1VarM, OtherVarM,
priorcode[1], DD))
if (nlev > 1 && !isTRUE(xc)) {
wrt("ERAS G30")
wrt("LINK 0 G30")
}
wrt("ERAS c1090 c1091")
wrt("")
if (!is.null(dami) && dami[1] == 0 && length(dami) > 1) {
ndami <- length(dami)
mvnames <- rep(NA, ndami - 1)
wrt(paste("LINK", ndami, "G23"))
for (i in 2:ndami) {
wrt(paste("MCMC 1 ", dami[i] - dami[i - 1], " ", thinning, " c1090 c1091 c1003 c1004 1 ", DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt(paste0("DAMI 0 G23[", i - 1, "]"))
mvnames[i - 1] <- paste0("'_est_", dami[i], "'")
wrt(paste0("NAME G23[", i - 1, "] ", mvnames[i - 1]))
wrt("PAUS 1")
}
wrt("LINK 0 G23")
if (dami[ndami] < iterations) {
wrt(paste("MCMC 1 ", (iterations/thinning) - dami[ndami], " ", thinning, " c1090 c1091 c1003 c1004 1 ",
DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt("PAUS 1")
}
wrt(paste("PSTA", paste0("'", MIfile, "'"), paste(mvnames, collapse = " "), "'resp_indicator'"))
wrt(paste("ERAS ", paste(mvnames, collapse = " "), sep = ""))
} else {
if (debugmode) {
for (i in 1:floor((iterations/thinning)/refresh)) {
wrt(paste("MCMC 1 ", refresh, " ", thinning, " c1090 c1091 c1003 c1004 1 ", DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt("PAUS 1")
}
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
if (!is.wholenumber(iterations/refresh)) {
wrt(paste("MCMC 1 ", (iterations/thinning)%%refresh, " ", thinning, " c1090 c1091 c1003 c1004 1 ",
DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt("PAUS 1")
}
} else {
wrt(paste("MCMC 1 ", iterations/thinning, " ", thinning, " c1090 c1091 c1003 c1004 1 ", DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
}
}
wrt("ECHO 0")
if (!is.null(saveworksheet)) {
wrt(paste0("STOR ", saveworksheet))
}
if (debugmode) {
wrt("WSET 15 1")
wrt("EXPA 3")
wrt("ESTM 2")
wrt("PAUS")
}
wrt("NOTE *****************************************************************")
wrt("NOTE Export the model results to R")
wrt("NOTE *****************************************************************")
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
wrt("NOBS 2 b31 b32")
} else {
wrt("NOBS 1 b31 b32")
}
wrt("EDIT 1 '_Stats' b31")
wrt("EDIT 2 '_Stats' b32")
if (!(D[1] == "Mixed") && is.null(merr) && is.null(fact)) {
wrt("BDIC b1 b2 b3 b4")
wrt("EDIT 3 '_Stats' b1")
wrt("EDIT 4 '_Stats' b2")
wrt("EDIT 5 '_Stats' b3")
wrt("EDIT 6 '_Stats' b4")
}
wrt("EDIT 7 '_Stats' b21")
wrt("EDIT 8 '_Stats' b22")
wrt("NAME c1098 '_FP_b'")
wrt("NAME c1099 '_FP_v'")
wrt("NAME c1096 '_RP_b'")
wrt("NAME c1097 '_RP_v'")
wrt("NAME c1094 '_esample'")
wrt("SUM '_esample' b1")
wrt("EDIT 9 '_Stats' b1")
if (is.null(fact)) {
wrt(paste("PSTA '", MCMCfile, "' ", "'_FP_b' ", "'_FP_v' ", "'_RP_b' ", "'_RP_v' ", "'_Stats'", sep = ""))
} else {
wrt("LINK 6 G25")
wrt("DAFA G25[1] G25[2]")
wrt("DAFL G25[3] G25[4]")
wrt("DAFV G25[5] G25[6]")
wrt("NAME G21[1] '_FACT_load_b_chain'")
wrt("NAME G21[2] '_FACT_load_v_chain'")
wrt("NAME G25[1] '_FACT_value_b'")
wrt("NAME G25[2] '_FACT_value_v'")
wrt("NAME G25[3] '_FACT_load_b'")
wrt("NAME G25[4] '_FACT_load_v'")
wrt("NAME G25[5] '_FACT_var_b'")
wrt("NAME G25[6] '_FACT_var_v'")
wrt(paste("PSTA '", FACTchainfile, "' ", "'_FACT_load_b_chain' ", "'_FACT_load_v_chain' ", "'_FACT_value_b' ",
"'_FACT_value_v' ", sep = ""))
wrt(paste("PSTA '", MCMCfile, "' ", "'_FP_b' ", "'_FP_v' ", "'_RP_b' ", "'_RP_v' ", "'_FACT_load_b' ",
"'_FACT_load_v' ", "'_FACT_var_b' ", "'_FACT_var_v' ", "'_Stats'", sep = ""))
wrt("ERAS G21")
wrt("LINK 0 G21")
wrt("ERAS G25")
wrt("LINK 0 G25")
}
wrt("ERAS '_Stats'")
wrt("")
if (!is.null(dami) && length(dami) == 1) {
wrt("NOTE save imputed values if there are missing values")
wrt("SWIT b1")
wrt("CASE 0:")
wrt("LEAVE")
wrt("CASE:")
wrt("LINK 3 G26")
wrt("NAME G26[1] '_MissingInd'")
wrt("CALC '_MissingInd'=abso('_esample'-1)")
if (dami == 1) {
wrt("DAMI 1 G26[2]")
wrt("NAME G26[2] '_est'")
wrt(paste("PSTA '", MIfile, "' '_est' '_MissingInd' ", sep = ""))
wrt("ERAS '_est'")
}
if (dami == 2) {
wrt("DAMI 2 G26[2] G26[3]")
wrt("NAME G26[2] '_est'")
wrt("NAME G26[3] '_SDs'")
wrt(paste("PSTA '", MIfile, "' '_est' '_SDs' '_MissingInd' ", sep = ""))
wrt("ERAS '_est' '_SDs'")
}
wrt("LINK 0 G26")
wrt("ENDS")
wrt("")
}
wrt("NOTE export parameter chain")
wrt("NAME c1091 'deviance'")
wrt("NAME c1090 'mcmcchains'")
wrt("LINK 0 G25")
wrt("LINK 0 G26")
if (D[1] == "Multinomial") {
nresp <- length(levels(indata[, resp])) - 1
resp.names <- levels(indata[, resp])[-as.numeric(D[5])]
if (is.list(expl)) {
nonfp.s <- nonfp.sep
for (i in 1:length(resp.names)) {
if (D["mode"] == 0) {
nonfp.s <- gsub(paste(".", resp.names[i], sep = ""), "", nonfp.s)
}
if (D["mode"] == 1) {
nonfp.s <- gsub(paste(".(>=", resp.names[i], ")", sep = ""), "", nonfp.s)
}
}
nonfp.s <- unique(nonfp.s)
if (is.na(sep.coeff[1]))
sep.coeff <- character(0)
for (p in sep.coeff) {
if (is.na(nonfp.sep[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
# svec.common <- character(0)
kk <- 1
tempid <- 1:(nresp + 1)
tempid <- tempid[-as.numeric(D["ref.cat"])]
for (p in common.coeff) {
newp <- paste(p, paste(tempid[as.logical(common.coeff.id[kk, ])], collapse = ""), sep = ".")
kk <- kk + 1
nonfp.c <- nonfp.common
if (is.na(nonfp.common[1]) || sum(newp == nonfp.c) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", newp)), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", newp), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
} else {
nonfp.s <- nonfp
for (i in 1:length(resp.names)) {
if (D["mode"] == 0) {
nonfp.s <- gsub(paste(".", resp.names[i], sep = ""), "", nonfp.s)
}
if (D["mode"] == 1) {
nonfp.s <- gsub(paste(".(>=", resp.names[i], ")", sep = ""), "", nonfp.s)
}
}
nonfp.s <- unique(nonfp.s)
expla <- expl
for (p in expla) {
if (is.na(nonfp[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
}
} else {
if (D[1] == "Multivariate Normal" || D[1] == "Mixed") {
nresp <- length(resp)
if (is.list(expl)) {
nonfp.s <- nonfp.sep
for (i in 1:length(resp)) {
nonfp.s <- gsub(paste(".", resp[i], sep = ""), "", nonfp.s)
}
nonfp.s <- unique(nonfp.s)
if (is.na(sep.coeff[1]))
sep.coeff <- character(0)
for (p in sep.coeff) {
if (is.na(nonfp.sep[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
kk <- 1
for (p in common.coeff) {
newp <- paste(p, paste(which(as.logical(common.coeff.id[kk, ])), collapse = ""), sep = ".")
kk <- kk + 1
nonfp.c <- nonfp.common
if (is.na(nonfp.common[1]) || sum(newp == nonfp.c) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", newp)), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", newp), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
} else {
nonfp.s <- nonfp
for (i in 1:length(resp)) {
nonfp.s <- gsub(paste(".", resp[i], sep = ""), "", nonfp.s)
}
nonfp.s <- unique(nonfp.s)
expla <- expl
for (p in expla) {
if (is.na(nonfp[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
}
} else {
expla <- expl
for (p in expla) {
if (is.na(nonfp[1]) || sum(p == nonfp) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", p), "'"))
wrt(paste0("DESC G25[1] 'FP:", p, "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
}
wrt.resid <- function(rpx, resid.lev) {
nrpx <- length(rpx)
for (j in 1:nrpx) {
for (i in 1:j) {
if (i == j) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_var_", chartr(".", "_", rpx[i])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":var(", chartr(".", "_", rpx[i]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_cov_", chartr(".", "_", rpx[i]), "_", chartr(".", "_", rpx[j])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":cov(", chartr(".", "_", rpx[i]), ",", chartr(".", "_", rpx[j]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
wrt.resid2 <- function(rpx, resid.lev, clre) {
nrpx <- length(rpx)
nclre <- ncol(clre)
k <- 1
for (j in 1:nrpx) {
for (i in 1:j) {
if (i == j) {
if (resid.lev == as.numeric(clre[1, k]) && rpx[i] == clre[2, k] && rpx[i] == clre[3, k]) {
if (k < ncol(clre))
k <- k + 1
} else {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_var_", chartr(".", "_", rpx[i])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":var(", chartr(".", "_", rpx[i]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if ((resid.lev == as.numeric(clre[1, k]) && rpx[i] == clre[2, k] && rpx[j] == clre[3, k]) || (resid.lev ==
as.numeric(clre[1, k]) && rpx[j] == clre[2, k] && rpx[i] == clre[3, k])) {
if (k < ncol(clre))
k <- k + 1
} else {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_cov_", chartr(".", "_", rpx[i]), "_", chartr(".", "_", rpx[j])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":cov(", chartr(".", "_", rpx[i]), ",", chartr(".", "_", rpx[j]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
wrt.resid3 <- function(rpx, resid.lev) {
nrpx <- length(rpx)
for (j in 1:nrpx) {
for (i in 1:j) {
if (i == j) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_var_", chartr(".", "_", rpx[i])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":var(", chartr(".", "_", rpx[i]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
if (nrp > 0) {
for (ii in 1:nrp) {
if (!is.null(fact)) {
wrt.resid3(rp[[ii]], as.numeric(sub("rp", "", rp.names[ii])))
} else {
if (is.null(clre)) {
wrt.resid(rp[[ii]], as.numeric(sub("rp", "", rp.names[ii])))
} else {
wrt.resid2(rp[[ii]], as.numeric(sub("rp", "", rp.names[ii])), clre)
}
}
}
}
# Add in extra parameters ect.
if (D[1] == "Multinomial" && as.numeric(D["mode"]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] 'RP1_bcons_1'"))
wrt(paste0("DESC G25[1] 'RP1:bcons_1'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
wrt("GSIZ G26 b1000")
wrt("LINK 2 G27")
wrt(paste("CODE b1000", 1, iterations/thinning, "G27[1]"))
wrt("NAME G27[1] 'parnum'")
wrt(paste("GENErate", 1, iterations/thinning, "G27[2]"))
wrt("NAME G27[2] 'iteration'")
wrt("DESC G27[2] '\\Iteration'")
wrt("SPLIt 'mcmcchains' 'parnum' G26")
if (D[1] == "Multinomial" && as.numeric(D["mode"]) == 1) {
wrt("LINK 1 G25")
wrt(paste("PUT ", iterations/thinning, 1, "G25[1]"))
wrt("NAME G25[1] 'RP1_bcons_1'")
wrt("DESC G25[1] 'RP1:bcons_1'")
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
wrt(paste0("PSTA '", chainfile, "' ", "'iteration' 'deviance' G26"))
wrt("ERAS 'parnum' 'iteration' G26")
wrt("LINK 0 G27")
wrt("LINK 0 G26")
calcresiduals <- function(level, displevel, rpx, resioptions, clre = clre) {
wrt("")
len.rpx <- length(rpx)
# For MCMC there is always only one residual column at level 1
if (level == 1) {
len.rpx = 1
}
ii <- 1
wrt(paste("LINK", len.rpx, "G26"))
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G26[", k, "] '", shortname("lev_", displevel, "_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G26[", k, "] ", "'residual estimates'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G26[", k, "] '", shortname("lev_", displevel, "_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G26[", k, "] ", "'residual estimates'"))
}
}
wrt(paste("LINK", len.rpx, "G27"))
if ("variance" %in% resioptions) {
ii <- 1
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_variance_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual variance'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_variance_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual variance'"))
}
}
} else {
residual_se <- NULL
ii <- 1
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_se_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual standard error'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_se_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual standard error'"))
}
}
}
wrt("RFUN")
wrt("ROUT G26 G27")
wrt("")
wrt(paste("RLEV ", level, sep = ""))
wrt("RCOV 1")
outgroups <- c("G26", "G27")
if ("standardised" %in% resioptions) {
wrt(paste("LINK", len.rpx, "G28"))
ii <- 1
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G28[", k, "] '", shortname("lev_", displevel, "_std_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G28[", k, "] ", "'std standardised residual'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G28[", k, "] '", shortname("lev_", displevel, "_std_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G28[", k, "] ", "'std standardised residual'"))
}
}
wrt("RTYP 0")
wrt("MCRE")
ccount <- 1
for (k in 1:len.rpx) {
wrt(paste0("CALC G28[", k, "]=G26[", k, "]/sqrt(G27[", k, "])", sep = ""))
}
outgroups <- c(outgroups, "G28")
}
wrt("RTYP 1") # Compute comparative variances
wrt("MCRE")
if (!("variance" %in% resioptions)) {
for (k in 1:len.rpx) {
wrt(paste0("CALC G27[", k, "]=sqrt(G27[", k, "])")) # Convert the variances to standard errors
}
}
wrt("")
wrt(paste0("NOBS ", level, " b30 b31"))
wrt("LINK 1 G29")
wrt("GENE 1 b30 1 G29[1]")
wrt(paste0("NAME G29[1] 'lev_", displevel, "_residualid'"))
outgroups <- c(outgroups, "G29")
wrt("")
wrt("LINK 0 G30")
for (i in 1:length(outgroups)) {
wrt(paste("GSET 2 G30", outgroups[i], "G30"))
wrt(paste("LINK 0", outgroups[i]))
}
}
if (resi.store && nrp > 0) {
for (j in nrp:1) {
rpx <- rp[[j]]
len.rpx <- length(rp[[j]])
wrt(paste("NOTE Calculate level ", as.numeric(sub("rp", "", rp.names[j])), " residuals", sep = ""))
levtt <- as.numeric(sub("rp", "", rp.names[j]))
displevel <- levtt
if (D[1] == "Multivariate Normal" || D[1] == "Mixed" || D[1] == "Multinomial") {
displevel <- displevel - 1
}
calcresiduals(levtt, displevel, rpx, resioptions, clre = clre)
wrt(paste("PSTA '", resifile[j], "' G30", sep = ""))
}
wrt("ERAS G30")
wrt("LINK 0 G30")
if (!is.null(resi.store.levs)) {
resiname <- rep(NA, length(resi.store.levs))
for (i in 1:length(resi.store.levs)) {
resiname[i] <- paste("'resi_lev", resi.store.levs[i], "'", sep = "")
}
wrt(paste("PSTA '", resichains, "' ", paste(resiname, collapse = " "), sep = ""))
}
}
}
if (debugmode) {
wrt("WPMT 'Do you want to close MLwiN?' b50")
wrt("SWITch b50")
wrt(" CASE 1:")
wrt(" EXIT ")
} else {
wrt("EXIT")
}
return(namemap)
}
Try the R2MLwiN package in your browser
Any scripts or data that you put into this service are public.
R2MLwiN documentation built on May 29, 2024, 2:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions write.MCMC
Documented in write.MCMC
#' Writes MLwiN macros to fit models using Markov chain Monte Carlo (MCMC)
#' methods
#'
#' write.MCMC is an internal function which creates an MLwiN macro file to
#' fit models using MCMC.
#'
#' @param indata A data.frame object containing the data to be modelled.
#' @param dtafile The file name of the dataset to be imported into MLwiN, which
#' is in Stata format (i.e. with extension .dta).
#' @param oldsyntax Specified as \code{FALSE} if new syntax has been used in
#' \code{Formula} object, else specified as \code{TRUE} (to enable
#' backcompatibility).
#' @param resp A character string (vector) of response variable(s).
#' @param 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.
#' \code{levID = c('level2', 'level1')} where \code{'level2'} is the higher
#' level).
#' @param expl A character string (vector) of explanatory (predictor)
#' variable(s).
#' @param rp A character string (vector) of random part of random variable(s).
#' @param D A character string/vector specifying the type of distribution to be modelled, which
#' can include \code{'Normal'} (the default), \code{'Binomial'}, \code{'Poisson'},
#' \code{'Unordered Multinomial'}, \code{'Ordered Multinomial'},
#' \code{'Multivariate Normal'}, or \code{'Mixed'}. In the case of the latter,
#' \code{'Mixed'} precedes the response types which also need to be listed in
#' \code{D}, e.g. \code{c('Mixed', 'Normal', 'Binomial')}; these need to be
#' be listed in the same order to which they are referred to in the
#' \code{Formula} object (see \code{\link{runMLwiN}}, \code{\link{Formula.translate}},
#' \code{\link{Formula.translate.compat}}. \code{'Mixed'} combinations can only consist of
#' \code{'Normal'} and \code{'Binomial'} for MCMC estimation.
#' @param 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).
#' \code{N = 0} specifies marginal quasi-likelihood
#' linearization (MQL), whilst \code{N = 1} specifies penalised quasi-
#' likelihood linearization (PQL); \code{M = 1} specifies first order
#' approximation, whilst \code{M = 2} specifies second order approximation.
#' \code{nonlinear = c(N = 0, M = 1)} by default. First order marginal
#' quasi-likelihood (MQL1) only option for single-level discrete response
#' models.
#' @param 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),
#' \code{NA}(s) if no reference group; the third row states the number of
#' categories for each variable.
#' @param notation Specifies the model subscript notation to be used in the
#' MLwiN equations window. \code{'class'} means no multiple subscripts, whereas
#' \code{'level'} has multiple subscripts.
#' @param nonfp Removes the fixed part of random variable(s). \code{NA} if no
#' variable is removed.
#' @param 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.
#' @param Meth Specifies the maximum likelihood estimation method to be used
#' when generating starting values via (R)IGLS.
#' If \code{Meth = 0} estimation method is set to RIGLS. If \code{Meth = 1}
#' estimation method is set to IGLS (the default setting). If \code{Meth} is
#' absent, alternate between IGLS and RIGLS.
#' @param merr A vector which sets-up measurement errors on predictor
#' variables. The first element \code{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 \code{demo(MCMCGuide14)} for an
#' example.
#' @param carcentre If CAR model (i.e. if \code{car} is non-\code{NULL}),
#' \code{carcentre = TRUE} mean-centres all random effects at that level.
#' @param 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 \code{startval = NULL}).
#' @param convtol When generating starting values via (R)IGLS, a numeric
#' value specifying the IGLS convergence criterion, as
#' specified in the \code{tol} option within \code{estoptions}, where
#' \code{startval = NULL}) (see \code{\link{runMLwiN}}). If value of
#' \code{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 \code{2} for m if not
#' otherwise specified.
#' @param seed An integer specifying the random seed in MLwiN.
#' @param iterations An integer specifying the number of iterations after
#' burn-in.
#' @param burnin An integer specifying length of the burn-in.
#' @param 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.
#' @param thinning An integer specifying the frequency with which successive
#' values in the Markov chain are stored. By default \code{thinning = 1}.
#' @param priorParam A vector specifying the informative priors used, as output
#' from \code{\link{prior2macro}}.
#' @param refresh An integer specifying how frequently the parameter estimates
#' are refreshed on the screen during iterations; only applies if
#' \code{debugmode = TRUE} in \code{estoptions}:
#' see \code{\link[R2MLwiN]{runMLwiN}}.
#' @param fixM Specifies the fixed effect method: \code{1} for Gibbs Sampling,
#' \code{2} for univariate MH Sampling and \code{3} for multivariate MH
#' Sampling.
#' @param residM Specifies the residual method: \code{1} for Gibbs Sampling,
#' \code{2} for univariate MH Sampling and \code{3} for multivariate MH
#' Sampling.
#' @param Lev1VarM Specifies the level 1 variance method: \code{1} for Gibbs
#' Sampling, \code{2} for univariate MH Sampling and \code{3} for multivariate
#' MH Sampling.
#' @param OtherVarM Specifies the variance method for other levels: \code{1}
#' for Gibbs Sampling, \code{2} for univariate MH Sampling and \code{3} for
#' multivariate MH Sampling.
#' @param adaption \code{adaption = 1} indicates adaptation is to be used;
#' \code{0} otherwise.
#' @param priorcode A vector indicating which default priors are to be used
#' for the variance parameters. It defaults to \code{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. \code{c(gamma = 1, shape = 0.5, scale = 0.2)}). Alternatively
#' \code{c(uniform = 1)} specifies Uniform priors on the variance scale. To allow
#' for back-compatibility with deprecated syntax used in versions of
#' \pkg{R2MLwiN} prior to 0.8-2, if \code{priorcode} is instead specified as
#' an integer, then \code{1} indicates that Gamma priors are used, whereas
#' \code{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.
#' @param rate An integer specifying the acceptance rate (as a percentage);
#' this command is ignored if \code{adaption = 0}.
#' @param tol An integer specifying tolerance (as a percentage) for the acceptance rate.
#' @param lclo This command toggles on/off the possible forms of complex level
#' 1 variation when using MCMC. \code{lclo = 0} expresses the level
#' 1 variation as a function of the predictors, whereas \code{lclo = 1} expresses the
#' log of the level 1 precision (1/variance) as a function of the predictors.
#' @param mcmcOptions A list of other MCMC options used. See `Details' below.
#' @param fact A list of objects specified for factor analysis. See `Details'
#' below.
#' @param xc Indicates whether model is cross-classified (\code{TRUE}) or
#' nested (\code{FALSE}). \code{xc = NULL} by default (corresponding to
#' \code{FALSE}), unless either \code{mm} or \code{car} are not null, in
#' which case \code{xc = TRUE}.
#' @param 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
#' \code{\link{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 \code{NA} is specified. Otherwise, lists need to be assigned to
#' \code{mmvar} and \code{weights}, with the former containing columns
#' specifying the classification units, and the latter containing columns
#' specifying the weights. Ignored if \code{EstM = 0}, i.e. only applicable to models estimated via
#' MCMC. \code{mm = NULL} by default. Supersedes deprecated \code{xclass}.
#' E.g. (from \code{demo(MCMCGuide16)}) for
#' \code{logearn ~ 1 + age_40 + sex + parttime + (company|1) + (id|1)}, if
#' \code{company} is a multiple membership classification with the variables
#' indicating the classifications in \code{company}, \code{company2},
#' \code{company3}, \code{company4} and their weights in \code{weight1}, \code{weight2},
#' \code{weight3} and \code{weight4} then
#' \code{mm = list(list(mmvar = list('company', 'company2', 'company3', 'company4'),}
#' \code{weights = list('weight1', 'weight2', 'weight3', 'weight4')), NA)}
#' with the \code{NA}, listed last, corresponding to the level 1 identifier (\code{id}).
#' @param 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 \code{NA} is specified. Otherwise, lists need to be assigned to
#' \code{carvar} and \code{weights}, with the former containing columns
#' specifying the spatial classification units, and the latter containing
#' columns specifying the weights. See \code{demo(MCMCGuide17)} for examples.
#' \code{car = NULL} by default.
#' @param BUGO If non-\code{NULL} uses BUGS for MCMC estimation using files
#' specified in \code{modelfile}, \code{initfile} and \code{datafile}.
#' @param 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 \code{nlev}+1; worksheet size in thousands of cells
#' is 6000; the number of columns is 2500; the number of explanatory variables
#' is \code{num_vars}+10; the number of group labels is 20. \code{nlev} is the
#' number of levels specified by \code{levID}, and \code{num_vars} is
#' approximately the number of explanatory variables calculated initially.
#' @param optimat This option instructs MLwiN to limit the maximum matrix size
#' that can be allocated by the (R)IGLS algorithm. Specify \code{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
#' \code{optimat = TRUE} caps the maximum matrix size at 800 lower-level units,
#' circumventing the MLwiN error message, and allowing most MLwiN
#' functionality.
#' @param modelfile A file name where the BUGS model will be saved in .txt
#' format.
#' @param initfile A file name where the BUGS initial values will be saved
#' in .txt format.
#' @param datafile A file name where the BUGS data will be saved in .txt
#' format.
#' @param macrofile A file name where the MLwiN macro file will be saved.
#' @param IGLSfile A file name where the IGLS estimates will be saved.
#' @param MCMCfile A file name where the MCMC estimates will be saved.
#' @param chainfile A file name where the MCMC chains will be saved.
#' @param MIfile A file name where the missing values will be saved.
#' @param resifile A file name where the residual estimates will be saved.
#' @param resi.store A logical value to indicate if residuals are to be stored
#' (\code{TRUE}) or not (\code{FALSE}).
#' @param resioptions A string vector to specify the various residual options.
#' The \code{'variance'} option calculates the posterior variances instead of
#' the posterior standard errors; the \code{'standardised'} option calculates standardised
#' residuals.
#' @param resichains A file name where the residual chains will be saved.
#' @param FACTchainfile A file name where the factor chains will be saved.
#' @param resi.store.levs An integer vector indicating the levels at which the
#' residual chains are to be stored.
#' @param debugmode A logical value determining whether MLwiN is run in the
#' background or not. The default value is \code{FALSE}: i.e. MLwiN is run in
#' the background. If \code{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 \code{workdir} ready to be read by \pkg{R2MLwiN}. Users can
#' instead opt to 'Abort macro' in which case the outputs are not saved to the
#' \code{workdir}. This option currently
#' works for 32 bit version of MLwiN only (automatically switches unless
#' \code{MLwiNPath} or \code{options(MLwiNPath)}
#' has been set directly to the executable).
#' @param startval A list of numeric vectors specifying the starting values
#' when using MCMC. \code{FP.b} corresponds to the estimates for the fixed
#' part; \code{FP.v} specifies the variance/covariance estimates for the fixed
#' part; \code{RP.b} specifies the variance estimates for the random part;
#' \code{RP.v} corresponds to the variance/covariance matrix of the variance
#' estimates for the random part.
#' @param dami This command outputs a complete (i.e. including non-missing
#' responses) response variable y. If \code{dami = c(0, <iter1>, <iter2>,...)} then
#' the response variables returned will be the value of y at the iterations
#' quoted (as integers \code{<iter1>, <iter2>}, etc.); these can be used for
#' multiple imputation. If \code{dami = 1} the value of y will be the mean
#' estimate from the iterations produced. \code{dami = 2} is as for \code{dami = 1}
#' but with the standard errors of the estimate additionally being stored.
#' @param namemap A mapping of column names to DTA friendly shorter names
#' @param saveworksheet A list of file names (one for each chain) used to store the
#' MLwiN worksheet after the model has been estimated.
#' @param rng.version An integer indicating which random number generator should
#' be used by MLwiN.
#'
#' @details A list of other MCMC options as used in the argument
#' \code{mcmcOptions}:
#' \itemize{
#' \item \code{orth}: If \code{orth = 1}, orthogonal fixed effect
#' vectors are used; zero otherwise.
#' \item \code{hcen}: An integer specifying the
#' level where we use hierarchical centering.
#' \item \code{smcm}: If \code{smcm = 1},
#' structured MCMC is used; zero otherwise.
#' \item \code{smvn}: If \code{smvn = 1}, the
#' structured MVN framework is used; zero otherwise.
#' \item \code{paex}: A matrix of Nx2; in each row, if the second digit is \code{1}, parameter expansion
#' is used at level <the first digit>.
#' \item \code{mcco}: This
#' command allows the user to have constrained settings for the lowest level
#' variance matrix in a multivariate Normal model. If value is \code{0},
#' it estimates distinct variances for each residual error and distinct covariances
#' for each residual error pair. Four other
#' settings are currently available:\cr
#' \tabular{ll}{\code{1} \tab fits stuctured errors with a common correlation paramater and a common variance parameter;\cr
#' \code{2} \tab fits AR1 errors with a common variance parameter;\cr \code{3} \tab fits structured errors with a common
#' correlation parameter and independent variance parameters;\cr \code{4} \tab fits AR1 errors with independent variance
#' parameters.\cr }
#' }
#'
#' A list of objects specified for cross-classified and/or multiple membership
#' models, as used in the argument \code{xclass}:
#' \itemize{
#' \item \code{class}: An integer
#' (vector) of the specified class(es).
#' \item \code{N1}: This defines a multiple
#' membership across \code{N1} units at level \code{class}. \code{N1}>1 if
#' there is multiple membership.
#' \item \code{weight}: If there is multiple
#' membership then the column number \code{weight}, which is the length of the
#' dataset, will contain the first set of weights for the multiple membership.
#' Note that there should be \code{N1} weight columns and they should be
#' sequential in the worksheet starting from \code{weight}.
#' \item \code{id}: If the
#' response is multivariate then the column number \code{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 \code{id} and following columns. There should be \code{N1} of
#' these identifier columns and they should be sequential starting from
#' \code{id} in the multivariate case.
#' \item \code{car}: \code{car = TRUE} indicates
#' the spatial CAR model; \code{FALSE} otherwise. \code{car = FALSE} if ignored.
#' }
#'
#' A list of objects specified for factor analysis, as used in the argument
#' \code{fact}:
#' \itemize{
#' \item \code{nfact}: Specifies the number of factors
#' \item \code{lev.fact}: Specifies the level/classification for the random part of
#' the factor for each factor.
#' \item \code{nfactcor}: Specifies the number of
#' correlated factors
#' \item \code{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
#' (\code{1}) or not (\code{0}). If more than one pair of factors is correlated,
#' then repeat this sequence for each pair.
#' \item \code{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.
#' \item \code{constr}: A matrix
#' specifying indicators of whether the factor loadings and the factor variance
#' are constrained (\code{1}) or not (\code{0}).
#' }
#'
#' @return Outputs a modified version of namemap containing newly generated
#' short names.
#' @note Note that for \code{FixM}, \code{residM}, \code{Lev1VarM} and
#'
#' \code{OtherVarM}, not all combinations of methods are available for all sets
#' of parameters and all models.
#'
#' @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.
#'
#' @author Zhang, Z., Charlton, C.M.J., Parker, R.M.A., Leckie, G., and Browne,
#' W.J. (2016) Centre for Multilevel Modelling, University of Bristol.
#'
#' @seealso \code{\link{write.IGLS}}
#'
write.MCMC <- function(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, rng.version = 10) {
shortname <- function(...) {
name <- paste0(...)
if (!name %in% names(namemap)) {
sname <- paste0("v", digest::digest(name, algo = "xxhash64", serialize = FALSE))
names(sname) <- name
namemap <<- c(namemap, sname)
}
return(namemap[[name]])
}
nlev <- length(levID)
nrp <- length(rp)
if (nrp > 0) {
rp.names <- names(rp)
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
for (i in 1:nrp) {
temp <- rp.names[i]
rp.names[i] <- paste("rp", as.integer(sub("rp", "", temp)) + 1, sep = "")
}
}
names(rp) <- rp.names
}
num.expl.init <- function(p, nonfp, categ) {
num_vars <- 0
if (is.na(nonfp[1])) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
num_vars <- num_vars + 1
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])])
} else {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1
}
}
} else {
if (sum(p == nonfp) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
num_vars <- num_vars + 1
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])])
} else {
num_vars <- num_vars + as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1
}
}
}
}
num_vars
}
if (is.list(expl)) {
sep.coeff <- expl$sep.coeff
if (is.list(nonfp)) {
nonfp.sep <- nonfp$nonfp.sep
nonfp.common <- nonfp$nonfp.common
}
if (!is.na(sep.coeff[1])) {
num_vars <- sum(sapply(sep.coeff, function(x) num.expl.init(x, nonfp.sep, categ)))
} else {
num_vars <- 0
}
if (D[1] == "Multinomial") {
nresp <- length(levels(indata[, resp])) - 1
## it is true if adding each expl variables separately
num_vars <- num_vars * nresp
}
if (D[1] == "Multivariate Normal") {
nresp <- length(resp)
num_vars <- num_vars * nresp
}
common.coeff <- expl$common.coeff
num_vars <- num_vars + sum(sapply(common.coeff, function(x) num.expl.init(x, nonfp$nonfp.common, categ)))
common.coeff.id <- expl$common.coeff.id
} else {
num_vars <- sum(sapply(expl, function(x) num.expl.init(x, nonfp, categ)))
if (D[1] == "Multinomial") {
nresp <- length(levels(indata[, resp])) - 1
num_vars <- num_vars * nresp
}
if (D[1] == "Multivariate Normal") {
nresp <- length(resp)
num_vars <- num_vars * nresp
}
}
if (nrp > 0) {
for (ii in 1:nrp) {
rpx <- rp[[rp.names[ii]]]
nrpx <- length(rpx)
if (nrpx == 1)
num_vars <- num_vars + 1
if (nrpx >= 2)
num_vars <- num_vars + nrpx * (nrpx - 1)/2 + nrpx
}
}
## Write into macro file
wrt <- function(x) write(x, macrofile, append = TRUE)
cat(file = macrofile)
wrt("ECHO 0")
wrt("NOTE *****************************************************************")
wrt("NOTE MLwiN macro created by rtomlwin command")
wrt("NOTE *****************************************************************")
wrt("")
wrt("NOTE Initialise MLwiN storage")
if (mem.init[1] == "default") {
wrt(paste("INIT ", nlev + 1, " 6000 2500 ", num_vars + 10, " 30", sep = ""))
} else {
wrt(paste("INIT ", mem.init[1], mem.init[2], mem.init[3], mem.init[4], mem.init[5]))
}
wrt("NOTE Limit the maximum matrix size")
if (optimat) {
wrt("OPTS 1")
} else {
wrt("OPTS 0")
}
wrt("MARK 0")
wrt("NOTE Import the R data into MLwiN")
wrt(paste("RSTA '", dtafile, "'", sep = ""))
wrt("NOTE Set the random number generator version")
wrt(paste0("RNGV ", rng.version))
wrt("NOTE Correct column names")
for (name in names(namemap)) {
wrt(paste0("COLN '", namemap[[name]], "' b50"))
wrt(paste0("DESC cb50 '", name, "'"))
wrt(paste0("NAME cb50 '", name, "'"))
}
if ("class" %in% notation) {
wrt("INDE 1")
}
if (!(D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed")) {
wrt("NOTE Specify the response variable")
for (p in resp) wrt(paste("RESP '", p, "'", sep = ""))
wrt("")
}
if (D[1] == "Binomial") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("RDISt 1 0")
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
if (D[2] == "logit") {
wrt("LFUN 0")
DD2 <- 0
}
if (D[2] == "probit") {
wrt("LFUN 1")
DD2 <- 1
}
if (D[2] == "cloglog") {
wrt("LFUN 2")
DD2 <- 2
}
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (D[1] == "Mixed") {
nresp <- length(resp)
for (ii in 1:nresp) wrt(paste("MVAR 1 '", resp[ii], "'", sep = ""))
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:2
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("IDEN 1 'resp_indicator'")
jj <- 1
for (ii in 2:length(D)) {
if (D[[ii]][1] == "Binomial") {
wrt(paste("RDISt ", jj, " 0", sep = ""))
wrt(paste("DOFFs ", jj, " '", D[[ii]][3], "'", sep = ""))
if (D[[ii]][2] == "logit") {
wrt("LFUN 0")
DD2 <- 0
}
if (D[[ii]][2] == "probit") {
wrt("LFUN 1")
DD2 <- 1
}
if (D[[ii]][2] == "cloglog") {
wrt("LFUN 2")
DD2 <- 2
}
}
if (D[[ii]][1] == "Poisson") {
wrt(paste("RDISt ", jj, " 1", sep = ""))
wrt("LFUN 3")
DD2 <- 3
DD2 <- 3
if (!is.na(D[[jj]][3])) {
wrt(paste("DOFFs ", jj, " '", D[[jj]][3], "'", sep = ""))
}
}
jj <- jj + 1
}
wrt("")
if (is.list(expl)) {
if (!is.na(sep.coeff[1])) {
interpos1 <- grep("\\:", sep.coeff)
if (!isTRUE(oldsyntax) || length(interpos1) == 0) {
for (x in 1:length(sep.coeff)) {
p <- sep.coeff[x]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- sep.coeff[interpos1]
explx <- sep.coeff[-interpos1]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
sep.coeff <- c(explx, exply)
}
}
interpos2 <- grep("\\:", common.coeff)
if (!isTRUE(oldsyntax) || length(interpos2) == 0) {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (!(y %in% interpos2)) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
} else {
TT <- ""
interx <- unlist(strsplit(p, "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
wrt("RPAT")
common.coeff[y] <- paste(p, ".", bb, sep = "")
}
}
} else {
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
}
wrt("")
}
if (D[1] == "Multivariate Normal") {
nresp <- length(resp)
for (ii in 1:nresp) wrt(paste("MVAR 1 '", resp[ii], "'", sep = ""))
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:2
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("IDEN 1 'resp_indicator'")
wrt("LFUN 0")
wrt("")
if (is.list(expl)) {
if (!is.na(sep.coeff[1])) {
interpos1 <- grep("\\:", sep.coeff)
if (!isTRUE(oldsyntax) || length(interpos1) == 0) {
for (x in 1:length(sep.coeff)) {
p <- sep.coeff[x]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- sep.coeff[interpos1]
explx <- sep.coeff[-interpos1]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
sep.coeff <- c(explx, exply)
}
}
interpos2 <- grep("\\:", common.coeff)
if (!isTRUE(oldsyntax) || length(interpos2) == 0) {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (!(y %in% interpos2)) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
} else {
TT <- ""
interx <- unlist(strsplit(p, "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
wrt("RPAT")
common.coeff[y] <- paste(p, ".", bb, sep = "")
}
}
} else {
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
}
wrt("")
}
if (D[1] == "Multinomial") {
wrt("LINK 2 G21")
wrt("NAME G21[1] 'resp' G21[2] 'resp_indicator'")
wrt("LINK 0 G21")
wrt(paste("MNOM ", as.numeric(D[4]), " '", resp, "' ", "'resp' 'resp_indicator' ", as.numeric(D[5]), sep = ""))
wrt("RESP 'resp'")
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:c(nlev - 1)) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("IDEN 1 'resp_indicator'")
wrt("")
if (as.numeric(D[4]) == 0)
wrt("RDISt 1 4") else wrt("RDISt 1 5")
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
if (D[2] == "logit") {
wrt("LFUN 0")
DD2 <- 0
}
if (D[2] == "probit") {
wrt("LFUN 1")
DD2 <- 1
}
if (D[2] == "cloglog") {
wrt("LFUN 2")
DD2 <- 2
}
wrt("")
if (is.list(expl)) {
if (!is.na(sep.coeff[1])) {
interpos1 <- grep("\\:", sep.coeff)
if (!isTRUE(oldsyntax) || length(interpos1) == 0) {
for (x in 1:length(sep.coeff)) {
p <- sep.coeff[x]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- sep.coeff[interpos1]
explx <- sep.coeff[-interpos1]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
sep.coeff <- c(explx, exply)
}
}
interpos2 <- grep("\\:", common.coeff)
if (!isTRUE(oldsyntax) || length(interpos2) == 0) {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
for (y in 1:length(common.coeff)) {
p <- common.coeff[y]
len.common.id <- length(common.coeff.id[y, ])
tt <- "RPAT "
aa <- 1:len.common.id
partname <- aa[rep(1, len.common.id) == common.coeff.id[y, ]]
bb <- ""
for (ii in aa) bb <- paste(bb, partname[ii], sep = "")
for (ii in 1:len.common.id) tt <- paste(tt, common.coeff.id[y, ii])
wrt(tt)
p <- common.coeff[y]
if (!(y %in% interpos2)) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
} else {
TT <- ""
interx <- unlist(strsplit(p, "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
wrt("RPAT")
common.coeff[y] <- paste(p, ".", bb, sep = "")
}
}
} else {
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
}
wrt("")
}
if (D[1] == "Normal") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("NOTE Specify covariate(s) used anywhere in the model")
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (D[1] == "Poisson") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("RDISt 1 1")
wrt("LFUN 3")
DD2 <- 3
if (!is.na(D[3])) {
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
}
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (D[1] == "Negbinom") {
wrt("NOTE Specify the level identifier(s)")
for (ii in 1:nlev) {
aa <- nlev:1
if (!is.na(levID[ii]))
wrt(paste("IDEN ", aa[ii], " '", levID[ii], "'", sep = ""))
}
wrt("")
wrt("RDISt 1 2")
wrt("LFUN 3")
DD2 <- 3
if (!is.na(D[3])) {
wrt(paste("DOFFs 1 '", D[3], "'", sep = ""))
}
interpos <- grep("\\:", expl)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in expl) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
} else {
exply <- expl[interpos]
explx <- expl[-interpos]
if (length(explx) > 0) {
for (p in explx) {
if (is.null(categ)) {
wrt(paste("ADDT '", p, "'", sep = ""))
} else {
if (sum(p == categ["var", ]) != 0) {
if (is.na(categ["ref", which(p == categ["var", ])])) {
wrt(paste("ADDT '", p, "' ", -1e+07, sep = ""))
} else {
wrt(paste("ADDT '", p, "' ", which(levels(indata[, p]) == categ["ref", which(p == categ["var",
])]), sep = ""))
}
} else {
wrt(paste("ADDT '", p, "'", sep = ""))
}
}
}
}
for (i in 1:length(exply)) {
TT <- ""
interx <- unlist(strsplit(exply[i], "\\:"))
for (j in 1:length(interx)) {
TT <- paste(TT, "'", interx[j], "' ", sep = "")
}
wrt(paste("ADDT ", TT, sep = ""))
}
expl <- c(explx, exply)
}
wrt("")
}
if (is.list(nonfp)) {
wrt("NOTE Turn off the fixed part of the explanatary variable(s)")
nonfp.sep <- nonfp$nonfp.sep
nonfp.common <- nonfp$nonfp.common
if (!is.na(nonfp.sep[1])) {
interpos <- grep("\\:", nonfp.sep)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in nonfp.sep) wrt(paste("FPAR 0 '", p, "'", sep = ""))
} else {
for (i in 1:length(nonfp.sep)) {
if (i %in% interpos && oldsyntax) {
wrt(paste("FPAR 0 '", gsub("\\:", "\\.", nonfp.sep[i]), "'", sep = ""))
} else {
wrt(paste("FPAR 0 '", nonfp.sep[i], "'", sep = ""))
}
}
}
}
if (!is.na(nonfp.common[1])) {
interpos <- grep("\\:", nonfp.common)
if (!isTRUE(oldsyntax) || length(interpos) == 0) {
for (p in nonfp.common) wrt(paste("FPAR 0 '", p, "'", sep = ""))
} else {
for (i in 1:length(nonfp.common)) {
if (i %in% interpos && oldsyntax) {
wrt(paste("FPAR 0 '", gsub("\\:", "\\.", nonfp.common[i]), "'", sep = ""))
} else {
wrt(paste("FPAR 0 '", nonfp.common[i], "'", sep = ""))
}
}
}
}
} else {
if (!is.na(nonfp[1])) {
wrt("NOTE Turn off the fixed part of the explanatory varible(s)")
for (p in nonfp) {
if(oldsyntax){
p <- gsub("\\:", "\\.", p)
}
wrt(paste("FPAR 0 '", p, "'", sep = ""))
}
}
}
wrt("")
wrt("NOTE Specify random part covariate(s)")
if (nrp > 0) {
for (ii in 1:nrp) {
for (p in rp[[ii]]) {
if (oldsyntax) {
p <- gsub("\\:", "\\.", p)
}
wrt(paste("SETV ", as.numeric(sub("rp", "", rp.names[ii])), " '", p, "'", sep = ""))
}
}
}
if (!is.null(clre)) {
nclre <- ncol(clre)
for (ii in 1:nclre) {
if (oldsyntax){
wrt(paste("CLRE ", as.numeric(clre[1, ii]), " '", gsub("\\:", "\\.", clre[2, ii]), "' '",
gsub("\\:", "\\.", clre[3, ii]), "'", sep = ""))
} else {
wrt(paste("CLRE ", as.numeric(clre[1, ii]), " '", clre[2, ii], "' '", clre[3, ii], "'", sep = ""))
}
}
}
nexpl <- length(expl)
wrt("")
wrt("NOTE Set estimation method")
if (Meth != 2) {
wrt(paste("METH", Meth))
}
wrt(paste("LINE ", nonlinear[1], nonlinear[2]))
wrt("")
if (!is.null(fact)) {
TT <- NULL
for (i in 1:fact$nfact) {
TT <- c(TT, fact$lev.fact[i] + 1, matrix(rbind(fact$loading[i, ], fact$constr[i, ]), nrow = 1))
}
if (fact$nfactcor > 0)
TT <- c(TT, fact$factcor)
FACT <- as.vector(c(length(resp), fact$nfact, fact$nfactcor, TT))
rm(TT)
TT <- ""
for (i in 1:length(FACT)) {
TT <- paste(TT, FACT[i])
}
wrt(paste("FACT ", TT))
wrt("LINK 2 G21")
wrt("SMFA 1 G21[1]")
wrt("SMFA 2 G21[2]")
}
if (D[1] == "Normal") {
wrt("PREF 0")
wrt("POST 0")
}
if ("simple" %in% notation) {
wrt("EXISt 'cons' b1000")
wrt("SWITch b1000")
wrt("CASE 0:")
wrt("EXISt 'Intercept' b1000")
wrt("SWITch b1000")
wrt("CASE 1:")
wrt("COLN 'Intercept' b1000")
wrt("NAME cb1000 'cons'")
wrt("ENDSWITch")
wrt("ENDSWITch")
wrt("NOTA 1")
}
wrt("NAME c1098 '_FP_b'")
wrt("NAME c1099 '_FP_v'")
wrt("NAME c1096 '_RP_b'")
wrt("NAME c1097 '_RP_v'")
wrt("NOTE Fit the model")
wrt("ECHO 1")
wrt("BATC 1")
wrt("MAXI 2")
wrt("STAR")
if (!is.null(startval)) {
if (!is.null(startval$FP.b)) {
wrt(paste("JOIN ", paste(startval$FP.b, collapse = " "), " '_FP_b'", sep = ""))
}
if (!is.null(startval$FP.v)) {
wrt(paste("JOIN ", paste(startval$FP.v[!upper.tri(startval$FP.v)], collapse = " "), " '_FP_v'", sep = ""))
}
if (!is.null(startval$RP.b)) {
wrt(paste("JOIN ", paste(startval$RP.b, collapse = " "), " '_RP_b'", sep = ""))
if (D[1] == "Multinomial" && D[4] == 0) {
wrt("JOIN '_RP_b' 1 '_RP_b'")
}
}
if (!is.null(startval$RP.v)) {
wrt(paste("JOIN ", paste(startval$RP.v[!upper.tri(startval$RP.v)], collapse = " "), " '_RP_v'", sep = ""))
if (D[1] == "Multinomial" && D[4] == 0) {
wrt(paste("JOIN '_RP_v' ", paste(rep(0, length(startval$RP.b)+1), collapse = " "), " '_RP_v'"))
}
}
} else {
wrt(paste("TOLE", convtol))
wrt(paste("MAXI", maxiter))
wrt("BATC 1")
wrt("NEXT")
}
wrt("ECHO 0")
wrt("ITNU 0 b21")
wrt("CONV b22")
wrt("")
wrt("NOTE *****************************************************************")
wrt("")
wrt("NOTE *****************************************************************")
wrt("NOTE Export the model results to R")
wrt("NOTE *****************************************************************")
wrt("LINK 1 G30")
wrt("NAME G30[1] '_Stats'")
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
wrt("NOBS 2 b31 b32")
} else {
wrt("NOBS 1 b31 b32")
}
wrt("EDIT 1 '_Stats' b31")
wrt("EDIT 2 '_Stats' b32")
if (D[1] == "Normal" || D[1] == "Multivariate Normal") {
wrt("LIKE b100")
}
wrt("EDIT 3 '_Stats' b100")
wrt("EDIT 7 '_Stats' b21")
wrt("EDIT 8 '_Stats' b22")
wrt("NAME c1094 '_esample'")
wrt("SUM '_esample' b1")
wrt("EDIT 9 '_Stats' b1")
wrt(paste("PSTA '", IGLSfile, "' ", "'_FP_b' ", "'_FP_v' ", "'_RP_b' ", "'_RP_v' ", "'_Stats'", sep = ""))
wrt("LINK 0 G30")
wrt("NOTE *****************************************************************")
wrt("NOTE Set estimation method to MCMC")
wrt("NOTE *****************************************************************")
wrt("EMODe 3")
if (!is.null(resi.store.levs)) {
wrt(paste0("LINK ", length(resi.store.levs), " G22"))
for (i in 1:length(resi.store.levs)) {
resilev <- resi.store.levs[i]
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
resilev <- resilev + 1
}
wrt(paste("SMRE ", resilev, " G22[", i, "]", sep = ""))
wrt(paste("NAME ", " G22[", i, "] 'resi_lev", resi.store.levs[i], "'", sep = ""))
}
}
if (!is.null(merr)) {
nmerr <- as.numeric(merr[1])
tt <- paste("MERR ", nmerr)
j <- 1
for (ii in 1:nmerr) {
tt <- paste(tt, " '", merr[j + 1], "' ", as.numeric(merr[j + 2]), sep = "")
j <- j + 2
}
wrt(tt)
}
wrt(paste("ORTH", mcmcOptions$orth))
if (mcmcOptions$hcen == 0) {
wrt(paste("HCEN", mcmcOptions$hcen))
} else {
wrt(paste("HCEN", 1, mcmcOptions$hcen))
}
wrt(paste("SMCM", mcmcOptions$smcm))
wrt(paste("SMVN", mcmcOptions$smvn))
if (is.matrix(mcmcOptions$paex)) {
apply(mcmcOptions$paex, 1, function(x) wrt(paste("PAEX", x[1], x[2])))
} else {
wrt(paste("PAEX", mcmcOptions$paex[1], mcmcOptions$paex[2]))
}
if (D[1] == "Multivariate Normal")
wrt(paste("MCCO ", mcmcOptions$mcco))
if (!is.null(xc) && isTRUE(xc)) {
wrt("XCLA 1")
}
if (!is.null(mm)) {
for (i in 1:length(mm)) {
if (!any(is.na(mm[[i]]))) {
if (D[1] != "Multivariate Normal") {
mmlev <- (length(mm) - i) + 1
wrt(paste0("MULM ", mmlev, " ", length(mm[[i]]$mmvar), " '", mm[[i]]$weights[[1]], "'"))
} else {
mmlev <- (length(mm) - i) + 2
wrt(paste0("MULM ", mmlev, " ", length(mm[[i]]$mmvar), " '", mm[[i]]$weights[[1]], "' '", mm[[i]]$mmvar[[1]],
"'"))
}
}
}
}
if (!is.null(car)) {
for (i in 1:length(car)) {
if (!any(is.na(car[[i]]))) {
carlev <- (length(car) - i) + 1
wrt(paste0("MULM ", carlev, " ", length(car[[i]]$carvar), " '", car[[i]]$weights[[1]], "' '", car[[i]]$carvar[[1]],
"'"))
wrt(paste("CARP", carlev, "1"))
}
}
if (!is.null(carcentre) && isTRUE(carcentre)) {
wrt("CARC 1")
}
}
wrt("")
wrt("NOTE Set MCMC seed")
wrt(paste("MCRS ", seed, sep = ""))
wrt("")
wrt("NOTE Set prior distribution parameters")
if (priorParam[1] != "default") {
wrt("PRIOR c1092")
lenpp <- length(priorParam)
tempt <- " "
for (i in 1:lenpp) tempt <- paste(tempt, priorParam[i])
wrt(paste("JOIN", tempt, " c1092", sep = ""))
wrt("")
} else {
wrt("")
}
if (D[1] == "Normal") {
wrt("PREF 0")
wrt("POST 0")
}
if (nlev > 1 && !isTRUE(xc)) {
wrt("MISR 0")
wrt("LINK 2 G30")
if (D[1] == "Multinomial" && as.numeric(D[4]) == 0) {
len.rpx <- 2
wrt(paste0("LINK ", len.rpx, " G27"))
wrt("LINK 1 G28")
wrt("NOTE Calculate MCMC starting values for level 2 residuals")
wrt("RLEV 2")
wrt("RFUN")
wrt("RCOV 2")
wrt("ROUT G27 G28")
wrt("RESI")
wrt("JOIN G30[1] G27 G30[1]")
wrt("JOIN G30[2] G28 G30[2]")
wrt("ERAS G27")
wrt("ERAS G28")
wrt("LINK 0 G27")
wrt("LINK 0 G28")
}
if (nrp > 0) {
for (j in nrp:1) {
if (as.numeric(sub("rp", "", rp.names[j])) != 1) {
rpx <- rp[[j]]
len.rpx <- length(rpx)
wrt(paste0("LINK ", len.rpx, " G27"))
wrt("LINK 1 G28")
wrt(paste("NOTE Calculate MCMC starting values for level ", as.numeric(sub("rp", "", rp.names[j])),
" residuals", sep = ""))
wrt(paste("RLEV ", as.numeric(sub("rp", "", rp.names[j])), sep = ""))
wrt("RFUN")
wrt("RCOV 2")
wrt("ROUT G27 G28")
wrt("RESI")
wrt("JOIN G30[1] G27 G30[1]")
wrt("JOIN G30[2] G28 G30[2]")
wrt("ERAS G27")
wrt("ERAS G28")
wrt("LINK 0 G27")
wrt("LINK 0 G28")
}
}
}
wrt("MISR 1")
}
if (D[1] == "Normal")
DD <- 1
if (D[1] == "Binomial")
DD <- 2
if (D[1] == "Mixed")
DD <- 5
if (D[1] == "Poisson")
DD <- 3
if (D[1] == "Multivariate Normal")
DD <- 4
if (D[1] == "Multinomial") {
if (as.numeric(D[4]) == 0)
DD <- 6 else DD <- 7
wrt("CLRV 2")
}
if (D[1] == "Negbinom")
DD <- 8
wrt(paste("LCLO ", lclo, sep = ""))
if (priorcode["gamma"] == 0) {
wrt("UNIP")
} else {
if (length(priorcode) > 1) {
wrt(paste("GAMP", priorcode["shape"], priorcode["scale"]))
}
}
if (debugmode) {
wrt("NOTE Open the equations window")
wrt("WSET 15 1")
wrt("EXPA 3")
wrt("ESTM 2")
wrt("PAUS")
}
priorcol <- ""
if (priorParam[1] != "default") {
priorcol <- "c1092"
}
if ((!is.null(BUGO)) && !(D[1] == "Mixed") && nrp > 0) {
if (D[1] == "Normal" || D[1] == "Multivariate Normal")
DD2 <- 0
if (nlev > 1 && !isTRUE(xc)) {
wrt(paste("BUGO 6 ", DD, " ", DD2, " G30[1] ", priorcol, " '", modelfile, "' ", "'", initfile, "' ", "'",
datafile, "'", sep = ""))
wrt("ERAS G30")
wrt("LINK 0 G30")
} else {
wrt(paste("BUGO 6 ", DD, " ", DD2, " ", "'", modelfile, "' ", "'", initfile, "' ", "'", datafile, "'",
sep = ""))
}
} else {
wrt("NOTE fit the model in MCMC")
wrt(paste("MTOT ", iterations, sep = ""))
wrt("ECHO 1")
residcols <- ""
if (nlev > 1 && !isTRUE(xc)) {
residcols <- "G30[1] G30[2]"
}
wrt(paste("MCMC 0", burnin, adaption, scale, rate, tol, residcols, priorcol, fixM, residM, Lev1VarM, OtherVarM,
priorcode[1], DD))
if (nlev > 1 && !isTRUE(xc)) {
wrt("ERAS G30")
wrt("LINK 0 G30")
}
wrt("ERAS c1090 c1091")
wrt("")
if (!is.null(dami) && dami[1] == 0 && length(dami) > 1) {
ndami <- length(dami)
mvnames <- rep(NA, ndami - 1)
wrt(paste("LINK", ndami, "G23"))
for (i in 2:ndami) {
wrt(paste("MCMC 1 ", dami[i] - dami[i - 1], " ", thinning, " c1090 c1091 c1003 c1004 1 ", DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt(paste0("DAMI 0 G23[", i - 1, "]"))
mvnames[i - 1] <- paste0("'_est_", dami[i], "'")
wrt(paste0("NAME G23[", i - 1, "] ", mvnames[i - 1]))
wrt("PAUS 1")
}
wrt("LINK 0 G23")
if (dami[ndami] < iterations) {
wrt(paste("MCMC 1 ", (iterations/thinning) - dami[ndami], " ", thinning, " c1090 c1091 c1003 c1004 1 ",
DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt("PAUS 1")
}
wrt(paste("PSTA", paste0("'", MIfile, "'"), paste(mvnames, collapse = " "), "'resp_indicator'"))
wrt(paste("ERAS ", paste(mvnames, collapse = " "), sep = ""))
} else {
if (debugmode) {
for (i in 1:floor((iterations/thinning)/refresh)) {
wrt(paste("MCMC 1 ", refresh, " ", thinning, " c1090 c1091 c1003 c1004 1 ", DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt("PAUS 1")
}
is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol
if (!is.wholenumber(iterations/refresh)) {
wrt(paste("MCMC 1 ", (iterations/thinning)%%refresh, " ", thinning, " c1090 c1091 c1003 c1004 1 ",
DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
wrt("PAUS 1")
}
} else {
wrt(paste("MCMC 1 ", iterations/thinning, " ", thinning, " c1090 c1091 c1003 c1004 1 ", DD, sep = ""))
wrt("PUPN c1003 c1004")
wrt("AVER c1091 b99 b100")
}
}
wrt("ECHO 0")
if (!is.null(saveworksheet)) {
wrt(paste0("STOR ", saveworksheet))
}
if (debugmode) {
wrt("WSET 15 1")
wrt("EXPA 3")
wrt("ESTM 2")
wrt("PAUS")
}
wrt("NOTE *****************************************************************")
wrt("NOTE Export the model results to R")
wrt("NOTE *****************************************************************")
if (D[1] == "Multinomial" || D[1] == "Multivariate Normal" || D[1] == "Mixed") {
wrt("NOBS 2 b31 b32")
} else {
wrt("NOBS 1 b31 b32")
}
wrt("EDIT 1 '_Stats' b31")
wrt("EDIT 2 '_Stats' b32")
if (!(D[1] == "Mixed") && is.null(merr) && is.null(fact)) {
wrt("BDIC b1 b2 b3 b4")
wrt("EDIT 3 '_Stats' b1")
wrt("EDIT 4 '_Stats' b2")
wrt("EDIT 5 '_Stats' b3")
wrt("EDIT 6 '_Stats' b4")
}
wrt("EDIT 7 '_Stats' b21")
wrt("EDIT 8 '_Stats' b22")
wrt("NAME c1098 '_FP_b'")
wrt("NAME c1099 '_FP_v'")
wrt("NAME c1096 '_RP_b'")
wrt("NAME c1097 '_RP_v'")
wrt("NAME c1094 '_esample'")
wrt("SUM '_esample' b1")
wrt("EDIT 9 '_Stats' b1")
if (is.null(fact)) {
wrt(paste("PSTA '", MCMCfile, "' ", "'_FP_b' ", "'_FP_v' ", "'_RP_b' ", "'_RP_v' ", "'_Stats'", sep = ""))
} else {
wrt("LINK 6 G25")
wrt("DAFA G25[1] G25[2]")
wrt("DAFL G25[3] G25[4]")
wrt("DAFV G25[5] G25[6]")
wrt("NAME G21[1] '_FACT_load_b_chain'")
wrt("NAME G21[2] '_FACT_load_v_chain'")
wrt("NAME G25[1] '_FACT_value_b'")
wrt("NAME G25[2] '_FACT_value_v'")
wrt("NAME G25[3] '_FACT_load_b'")
wrt("NAME G25[4] '_FACT_load_v'")
wrt("NAME G25[5] '_FACT_var_b'")
wrt("NAME G25[6] '_FACT_var_v'")
wrt(paste("PSTA '", FACTchainfile, "' ", "'_FACT_load_b_chain' ", "'_FACT_load_v_chain' ", "'_FACT_value_b' ",
"'_FACT_value_v' ", sep = ""))
wrt(paste("PSTA '", MCMCfile, "' ", "'_FP_b' ", "'_FP_v' ", "'_RP_b' ", "'_RP_v' ", "'_FACT_load_b' ",
"'_FACT_load_v' ", "'_FACT_var_b' ", "'_FACT_var_v' ", "'_Stats'", sep = ""))
wrt("ERAS G21")
wrt("LINK 0 G21")
wrt("ERAS G25")
wrt("LINK 0 G25")
}
wrt("ERAS '_Stats'")
wrt("")
if (!is.null(dami) && length(dami) == 1) {
wrt("NOTE save imputed values if there are missing values")
wrt("SWIT b1")
wrt("CASE 0:")
wrt("LEAVE")
wrt("CASE:")
wrt("LINK 3 G26")
wrt("NAME G26[1] '_MissingInd'")
wrt("CALC '_MissingInd'=abso('_esample'-1)")
if (dami == 1) {
wrt("DAMI 1 G26[2]")
wrt("NAME G26[2] '_est'")
wrt(paste("PSTA '", MIfile, "' '_est' '_MissingInd' ", sep = ""))
wrt("ERAS '_est'")
}
if (dami == 2) {
wrt("DAMI 2 G26[2] G26[3]")
wrt("NAME G26[2] '_est'")
wrt("NAME G26[3] '_SDs'")
wrt(paste("PSTA '", MIfile, "' '_est' '_SDs' '_MissingInd' ", sep = ""))
wrt("ERAS '_est' '_SDs'")
}
wrt("LINK 0 G26")
wrt("ENDS")
wrt("")
}
wrt("NOTE export parameter chain")
wrt("NAME c1091 'deviance'")
wrt("NAME c1090 'mcmcchains'")
wrt("LINK 0 G25")
wrt("LINK 0 G26")
if (D[1] == "Multinomial") {
nresp <- length(levels(indata[, resp])) - 1
resp.names <- levels(indata[, resp])[-as.numeric(D[5])]
if (is.list(expl)) {
nonfp.s <- nonfp.sep
for (i in 1:length(resp.names)) {
if (D["mode"] == 0) {
nonfp.s <- gsub(paste(".", resp.names[i], sep = ""), "", nonfp.s)
}
if (D["mode"] == 1) {
nonfp.s <- gsub(paste(".(>=", resp.names[i], ")", sep = ""), "", nonfp.s)
}
}
nonfp.s <- unique(nonfp.s)
if (is.na(sep.coeff[1]))
sep.coeff <- character(0)
for (p in sep.coeff) {
if (is.na(nonfp.sep[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
# svec.common <- character(0)
kk <- 1
tempid <- 1:(nresp + 1)
tempid <- tempid[-as.numeric(D["ref.cat"])]
for (p in common.coeff) {
newp <- paste(p, paste(tempid[as.logical(common.coeff.id[kk, ])], collapse = ""), sep = ".")
kk <- kk + 1
nonfp.c <- nonfp.common
if (is.na(nonfp.common[1]) || sum(newp == nonfp.c) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", newp)), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", newp), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
} else {
nonfp.s <- nonfp
for (i in 1:length(resp.names)) {
if (D["mode"] == 0) {
nonfp.s <- gsub(paste(".", resp.names[i], sep = ""), "", nonfp.s)
}
if (D["mode"] == 1) {
nonfp.s <- gsub(paste(".(>=", resp.names[i], ")", sep = ""), "", nonfp.s)
}
}
nonfp.s <- unique(nonfp.s)
expla <- expl
for (p in expla) {
if (is.na(nonfp[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
}
} else {
if (D[1] == "Multivariate Normal" || D[1] == "Mixed") {
nresp <- length(resp)
if (is.list(expl)) {
nonfp.s <- nonfp.sep
for (i in 1:length(resp)) {
nonfp.s <- gsub(paste(".", resp[i], sep = ""), "", nonfp.s)
}
nonfp.s <- unique(nonfp.s)
if (is.na(sep.coeff[1]))
sep.coeff <- character(0)
for (p in sep.coeff) {
if (is.na(nonfp.sep[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
kk <- 1
for (p in common.coeff) {
newp <- paste(p, paste(which(as.logical(common.coeff.id[kk, ])), collapse = ""), sep = ".")
kk <- kk + 1
nonfp.c <- nonfp.common
if (is.na(nonfp.common[1]) || sum(newp == nonfp.c) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", newp)), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", newp), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
} else {
nonfp.s <- nonfp
for (i in 1:length(resp)) {
nonfp.s <- gsub(paste(".", resp[i], sep = ""), "", nonfp.s)
}
nonfp.s <- unique(nonfp.s)
expla <- expl
for (p in expla) {
if (is.na(nonfp[1]) || sum(p == nonfp.s) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
for (j in 1:nresp) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", p), "_", resp[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", p), "_", resp[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (j in 1:nresp) {
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (j in 1:nresp) {
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i]), "_", resp.names[j]), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "_", resp.names[j], "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
}
} else {
expla <- expl
for (p in expla) {
if (is.na(nonfp[1]) || sum(p == nonfp) == 0) {
if (is.null(categ) || sum(p == categ["var", ]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", p), "'"))
wrt(paste0("DESC G25[1] 'FP:", p, "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
if (is.na(categ["ref", which(p == categ["var", ])])) {
categ.names <- levels(indata[[p]])
for (i in 1:as.numeric(categ["ncateg", which(p == categ["var", ])])) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
categ.names <- levels(indata[[p]])
refx <- categ["ref", which(p == categ["var", ])]
categ.names <- categ.names[-which(refx == categ.names)]
for (i in 1:(as.numeric(categ["ncateg", which(p == categ["var", ])]) - 1)) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("FP_", chartr(".", "_", categ.names[i])), "'"))
wrt(paste0("DESC G25[1] 'FP:", chartr(".", "_", categ.names[i]), "'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
}
}
wrt.resid <- function(rpx, resid.lev) {
nrpx <- length(rpx)
for (j in 1:nrpx) {
for (i in 1:j) {
if (i == j) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_var_", chartr(".", "_", rpx[i])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":var(", chartr(".", "_", rpx[i]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
} else {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_cov_", chartr(".", "_", rpx[i]), "_", chartr(".", "_", rpx[j])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":cov(", chartr(".", "_", rpx[i]), ",", chartr(".", "_", rpx[j]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
wrt.resid2 <- function(rpx, resid.lev, clre) {
nrpx <- length(rpx)
nclre <- ncol(clre)
k <- 1
for (j in 1:nrpx) {
for (i in 1:j) {
if (i == j) {
if (resid.lev == as.numeric(clre[1, k]) && rpx[i] == clre[2, k] && rpx[i] == clre[3, k]) {
if (k < ncol(clre))
k <- k + 1
} else {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_var_", chartr(".", "_", rpx[i])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":var(", chartr(".", "_", rpx[i]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
} else {
if ((resid.lev == as.numeric(clre[1, k]) && rpx[i] == clre[2, k] && rpx[j] == clre[3, k]) || (resid.lev ==
as.numeric(clre[1, k]) && rpx[j] == clre[2, k] && rpx[i] == clre[3, k])) {
if (k < ncol(clre))
k <- k + 1
} else {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_cov_", chartr(".", "_", rpx[i]), "_", chartr(".", "_", rpx[j])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":cov(", chartr(".", "_", rpx[i]), ",", chartr(".", "_", rpx[j]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
}
wrt.resid3 <- function(rpx, resid.lev) {
nrpx <- length(rpx)
for (j in 1:nrpx) {
for (i in 1:j) {
if (i == j) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] '", shortname("RP", resid.lev, "_var_", chartr(".", "_", rpx[i])), "'"))
wrt(paste0("DESC G25[1] 'RP", resid.lev, ":var(", chartr(".", "_", rpx[i]), ")'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
}
}
}
if (nrp > 0) {
for (ii in 1:nrp) {
if (!is.null(fact)) {
wrt.resid3(rp[[ii]], as.numeric(sub("rp", "", rp.names[ii])))
} else {
if (is.null(clre)) {
wrt.resid(rp[[ii]], as.numeric(sub("rp", "", rp.names[ii])))
} else {
wrt.resid2(rp[[ii]], as.numeric(sub("rp", "", rp.names[ii])), clre)
}
}
}
}
# Add in extra parameters ect.
if (D[1] == "Multinomial" && as.numeric(D["mode"]) == 0) {
wrt("LINK 1 G25")
wrt(paste0("NAME G25[1] 'RP1_bcons_1'"))
wrt(paste0("DESC G25[1] 'RP1:bcons_1'"))
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
wrt("GSIZ G26 b1000")
wrt("LINK 2 G27")
wrt(paste("CODE b1000", 1, iterations/thinning, "G27[1]"))
wrt("NAME G27[1] 'parnum'")
wrt(paste("GENErate", 1, iterations/thinning, "G27[2]"))
wrt("NAME G27[2] 'iteration'")
wrt("DESC G27[2] '\\Iteration'")
wrt("SPLIt 'mcmcchains' 'parnum' G26")
if (D[1] == "Multinomial" && as.numeric(D["mode"]) == 1) {
wrt("LINK 1 G25")
wrt(paste("PUT ", iterations/thinning, 1, "G25[1]"))
wrt("NAME G25[1] 'RP1_bcons_1'")
wrt("DESC G25[1] 'RP1:bcons_1'")
wrt("GSET 2 G26 G25 G26")
wrt("LINK 0 G25")
}
wrt(paste0("PSTA '", chainfile, "' ", "'iteration' 'deviance' G26"))
wrt("ERAS 'parnum' 'iteration' G26")
wrt("LINK 0 G27")
wrt("LINK 0 G26")
calcresiduals <- function(level, displevel, rpx, resioptions, clre = clre) {
wrt("")
len.rpx <- length(rpx)
# For MCMC there is always only one residual column at level 1
if (level == 1) {
len.rpx = 1
}
ii <- 1
wrt(paste("LINK", len.rpx, "G26"))
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G26[", k, "] '", shortname("lev_", displevel, "_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G26[", k, "] ", "'residual estimates'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G26[", k, "] '", shortname("lev_", displevel, "_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G26[", k, "] ", "'residual estimates'"))
}
}
wrt(paste("LINK", len.rpx, "G27"))
if ("variance" %in% resioptions) {
ii <- 1
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_variance_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual variance'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_variance_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual variance'"))
}
}
} else {
residual_se <- NULL
ii <- 1
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_se_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual standard error'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G27[", k, "] '", shortname("lev_", displevel, "_resi_se_", rpx[k]), "'"))
wrt(paste0("DESC G27[", k, "] ", "'residual standard error'"))
}
}
}
wrt("RFUN")
wrt("ROUT G26 G27")
wrt("")
wrt(paste("RLEV ", level, sep = ""))
wrt("RCOV 1")
outgroups <- c("G26", "G27")
if ("standardised" %in% resioptions) {
wrt(paste("LINK", len.rpx, "G28"))
ii <- 1
for (k in 1:len.rpx) {
if (!is.null(clre)) {
if (!(level == as.numeric(clre[1, ii]) && rpx[k] == clre[2, ii] && rpx[k] == clre[3, ii])) {
wrt(paste0("NAME G28[", k, "] '", shortname("lev_", displevel, "_std_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G28[", k, "] ", "'std standardised residual'"))
} else {
if (ii < ncol(clre))
ii <- ii + 1
}
} else {
wrt(paste0("NAME G28[", k, "] '", shortname("lev_", displevel, "_std_resi_est_", rpx[k]), "'"))
wrt(paste0("DESC G28[", k, "] ", "'std standardised residual'"))
}
}
wrt("RTYP 0")
wrt("MCRE")
ccount <- 1
for (k in 1:len.rpx) {
wrt(paste0("CALC G28[", k, "]=G26[", k, "]/sqrt(G27[", k, "])", sep = ""))
}
outgroups <- c(outgroups, "G28")
}
wrt("RTYP 1") # Compute comparative variances
wrt("MCRE")
if (!("variance" %in% resioptions)) {
for (k in 1:len.rpx) {
wrt(paste0("CALC G27[", k, "]=sqrt(G27[", k, "])")) # Convert the variances to standard errors
}
}
wrt("")
wrt(paste0("NOBS ", level, " b30 b31"))
wrt("LINK 1 G29")
wrt("GENE 1 b30 1 G29[1]")
wrt(paste0("NAME G29[1] 'lev_", displevel, "_residualid'"))
outgroups <- c(outgroups, "G29")
wrt("")
wrt("LINK 0 G30")
for (i in 1:length(outgroups)) {
wrt(paste("GSET 2 G30", outgroups[i], "G30"))
wrt(paste("LINK 0", outgroups[i]))
}
}
if (resi.store && nrp > 0) {
for (j in nrp:1) {
rpx <- rp[[j]]
len.rpx <- length(rp[[j]])
wrt(paste("NOTE Calculate level ", as.numeric(sub("rp", "", rp.names[j])), " residuals", sep = ""))
levtt <- as.numeric(sub("rp", "", rp.names[j]))
displevel <- levtt
if (D[1] == "Multivariate Normal" || D[1] == "Mixed" || D[1] == "Multinomial") {
displevel <- displevel - 1
}
calcresiduals(levtt, displevel, rpx, resioptions, clre = clre)
wrt(paste("PSTA '", resifile[j], "' G30", sep = ""))
}
wrt("ERAS G30")
wrt("LINK 0 G30")
if (!is.null(resi.store.levs)) {
resiname <- rep(NA, length(resi.store.levs))
for (i in 1:length(resi.store.levs)) {
resiname[i] <- paste("'resi_lev", resi.store.levs[i], "'", sep = "")
}
wrt(paste("PSTA '", resichains, "' ", paste(resiname, collapse = " "), sep = ""))
}
}
}
if (debugmode) {
wrt("WPMT 'Do you want to close MLwiN?' b50")
wrt("SWITch b50")
wrt(" CASE 1:")
wrt(" EXIT ")
} else {
wrt("EXIT")
}
return(namemap)
}
Try the R2MLwiN package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.