Nothing
#' @title Main function for GEM models
#' @description Main function to estimate and validate GEneralized Mixture models with uncertainty.
#' @aliases GEM
#' @usage GEM(Formula,family=c("cub","cube","ihg","cush"),data,...)
#' @param Formula Object of class Formula. Response variable is the vector of ordinal observations - see Details.
#' @param family Character string indicating which class of GEM models to fit.
#' @param data an optional data frame (or object coercible by \code{as.data.frame} to a data frame)
#' containing the variables in the model. If missing, the variables are taken from \code{environment(Formula)}.
#' @param ... Additional arguments to be passed for the specification of the model. See details and examples.
#' @export GEM
#' @return An object of the class "GEM" is a list containing the following elements:
#' \item{estimates}{Maximum likelihood estimates of parameters}
#' \item{loglik}{Log-likelihood function at the final estimates}
#' \item{varmat}{Variance-covariance matrix of final estimates}
#' \item{niter}{Number of executed iterations}
#' \item{BIC}{BIC index for the estimated model}
#' \item{ordinal}{Vector of ordinal responses on which the model has been fitted}
#' \item{time}{Processor time for execution}
#' \item{ellipsis}{Retrieve the arguments passed to the call and extra arguments generated via the call}
#' \item{family}{Character string indicating the sub-class of the fitted model}
#' \item{formula}{Returns the Formula of the call for the fitted model}
#' \item{call}{Returns the executed call}
#' @import methods Formula
#' @details It is the main function for GEM models estimation, calling for the corresponding function for
#' the specified subclass. The number of categories \code{m} is internally retrieved but it is advisable to pass
#' it as an argument to the call if some category has zero frequency.\cr
#' If \code{family="cub"}, then a CUB mixture model is fitted to the data to explain uncertainty,
#' feeling and possible shelter effect by further passing the extra argument \code{shelter} for the corresponding category.
#' Subjects' covariates can be included by specifying covariates matrices in the
#' Formula as \code{ordinal~Y|W|X}, to explain uncertainty (Y), feeling (W) or shelter (X). Notice that
#' covariates for shelter effect can be included only if specified for both feeling and uncertaint (GeCUB models). \cr
#' If \code{family="cube"}, then a CUBE mixture model (Combination of Uniform and Beta-Binomial) is fitted to the data
#' to explain uncertainty, feeling and overdispersion. Subjects' covariates can be also included to explain the
#' feeling component or all the three components by specifying covariates matrices in the Formula as
#' \code{ordinal~Y|W|Z} to explain uncertainty (Y), feeling (W) or
#' overdispersion (Z). An extra logical argument \code{expinform} indicates whether or not to use the expected or the
#' observed information matrix (default is FALSE). \cr
#' If \code{family="ihg"}, then an IHG model is fitted to the data. IHG models (Inverse Hypergeometric) are nested into
#' CUBE models (see the references below). The parameter \eqn{\theta} gives the probability of observing
#' the first category and is therefore a direct measure of preference, attraction, pleasantness toward the
#' investigated item. This is the reason why \eqn{\theta} is customarily referred to as the preference parameter of the
#' IHG model. Covariates for the preference parameter \eqn{\theta} have to be specified in matrix form in the Formula as
#' \code{ordinal~U}. \cr
#' If \code{family="cush"}, then a CUSH model is fitted to the data (Combination of Uniform and SHelter effect).
#' The category corresponding to the inflation should be
#' passed via argument \code{shelter}. Covariates for the shelter parameter \eqn{\delta}
#' are specified in matrix form Formula as \code{ordinal~X}. \cr
#' Even if no covariate is included in the model for a given component, the corresponding model matrix needs always
#' to be specified: in this case, it should be set to 0 (see examples below). Extra arguments include the maximum
#' number of iterations (\code{maxiter}, default: \code{maxiter}=500) for the optimization algorithm and
#' the required error tolerance (\code{toler}, default: \code{toler}=1e-6). \cr
#' Standard methods: \code{logLik()}, \code{BIC()}, \code{vcov()}, \code{fitted()}, \code{coef()}, \code{print()}, \code{summary()}
#' are implemented.\cr
#' The optimization procedure is run via \code{optim()} when required. If the estimated variance-covariance matrix is not positive definite, the function returns a
#' warning message and produces a matrix with NA entries.
#' @references
#' D'Elia A. (2003). Modelling ranks using the inverse hypergeometric distribution,
#' \emph{Statistical Modelling: an International Journal}, \bold{3}, 65--78 \cr
#' D'Elia A. and Piccolo D. (2005). A mixture model for preferences data analysis,
#' \emph{Computational Statistics & Data Analysis}, \bold{49}, 917--937 \cr
#' Capecchi S. and Piccolo D. (2017). Dealing with heterogeneity in ordinal responses,
#' \emph{Quality and Quantity}, \bold{51}(5), 2375--2393 \cr
#' Iannario M. (2014). Modelling Uncertainty and Overdispersion in Ordinal Data,
#' \emph{Communications in Statistics - Theory and Methods}, \bold{43}, 771--786 \cr
#' Piccolo D. (2015). Inferential issues for CUBE models with covariates,
#' \emph{Communications in Statistics. Theory and Methods}, \bold{44}(23), 771--786. \cr
#' Iannario M. (2015). Detecting latent components in ordinal data with overdispersion by means
#' of a mixture distribution, \emph{Quality & Quantity}, \bold{49}, 977--987 \cr
#' Iannario M. and Piccolo D. (2016a). A comprehensive framework for regression models of ordinal data.
#' \emph{Metron}, \bold{74}(2), 233--252.\cr
#' Iannario M. and Piccolo D. (2016b). A generalized framework for modelling ordinal data.
#' \emph{Statistical Methods and Applications}, \bold{25}, 163--189.\cr
#' @seealso \code{\link{logLik}}, \code{\link{coef}}, \code{\link{BIC}}, \code{\link{makeplot}},
#' \code{\link{summary}}, \code{\link{vcov}}, \code{\link{fitted}}, \code{\link{cormat}}
#' @keywords models
#' @examples
#'
#' library(CUB)
#' ## CUB models with no covariates
#' model<-GEM(Formula(Walking~0|0|0),family="cub",data=relgoods)
#' coef(model,digits=5) # Estimated parameter vector (pai,csi)
#' logLik(model) # Log-likelihood function at ML estimates
#' vcov(model,digits=4) # Estimated Variance-Covariance matrix
#' cormat(model) # Parameter Correlation matrix
#' fitted(model) # Fitted probability distribution
#' makeplot(model)
#' ################
#' ## CUB model with shelter effect
#' model<-GEM(Formula(officeho~0|0|0),family="cub",shelter=7,data=univer)
#' BICshe<-BIC(model,digits=4)
#' ################
#' ## CUB model with covariate for uncertainty
#' modelcovpai<-GEM(Formula(Parents~Smoking|0|0),family="cub",data=relgoods)
#' fitted(modelcovpai)
#' makeplot(modelcovpai)
#' ################
#' ## CUB model with covariates for both uncertainty and feeling components
#' data(univer)
#' model<-GEM(Formula(global~gender|freqserv|0),family="cub",data=univer,maxiter=50,toler=1e-2)
#' param<-coef(model)
#' bet<-param[1:2] # ML estimates of coefficients for uncertainty covariate: gender
#' gama<-param[3:4] # ML estimates of coefficients for feeling covariate: lage
#' ##################
#' ## CUBE models with no covariates
#' model<-GEM(Formula(MeetRelatives~0|0|0),family="cube",starting=c(0.5,0.5,0.1),
#' data=relgoods,expinform=TRUE,maxiter=50,toler=1e-2)
#' coef(model,digits=4) # Final ML estimates
#' vcov(model)
#' fitted(model)
#' makeplot(model)
#' summary(model)
#' ##################
#' ## IHG with covariates
#' modelcov<-GEM(willingn~freqserv,family="ihg",data=univer)
#' omega<-coef(modelcov) ## ML estimates
#' maxlik<-logLik(modelcov) ##
#' makeplot(modelcov)
#' summary(modelcov)
#' ###################
#' ## CUSH models without covariate
#' model<-GEM(Dog~0,family="cush",shelter=1,data=relgoods)
#' delta<-coef(model) # ML estimates of delta
#' maxlik<-logLik(model) # Log-likelihood at ML estimates
#' summary(model)
#' makeplot(model)
GEM<-function(Formula,family=c("cub","cube","ihg","cush"),data,...){
call <- match.call()
if (missing(data)) data <- environment(Formula)
mf<-model.frame(Formula,data,na.action=na.omit)
ordinal<-as.numeric(model.response(mf))
#if (missing(m)) m <- length(levels(factor(ordinal,ordered=TRUE)))
ellipsis.arg<-list(...)
ellipsis.arg[['data']]<-data
if (is.null(ellipsis.arg[['m']])){
ellipsis.arg[['m']] <- length(levels(factor(ordinal,ordered=TRUE)))
}
# m<-mget('m',ifnotfound=list(max(ordinal)))
### nc number of categories
#as.list(mapply(c, lista1, m, SIMPLIFY = TRUE))
# # #cl<-match.call(expand.dots=TRUE)
if (family == "cub"){
if (is.null(ellipsis.arg[['maxiter']])){
ellipsis.arg[['maxiter']] <- 500
}
if (is.null(ellipsis.arg$toler)){
ellipsis.arg$toler <- 1e-4
}
mod<-CUB(Formula,data,ellipsis.arg)
modelname<-"CUB"
}
if (family =="cube"){
if (is.null(ellipsis.arg[['maxiter']])){
ellipsis.arg[['maxiter']] <- 1000
}
if (is.null(ellipsis.arg$toler)){
ellipsis.arg$toler <- 1e-6
}
maxiter<- ellipsis.arg[['maxiter']]
toler<-ellipsis.arg$toler
mod<-CUBE(Formula,data,ellipsis.arg)
modelname<-"CUBE"
}
if (family == "ihg"){
# mod<-IHG(Formula,data,...)
modelname<-"IHG"
ellipsis.arg[['maxiter']]<-1
ellipsis.arg$niter<-1
mod<-IHG(Formula,data,ellipsis.arg)
}
if (family == "cush"){
ellipsis.arg[['maxiter']]<-1
ellipsis.arg$niter<-1
mod<-CUSH(Formula,data,ellipsis.arg)
modelname<-"CUSH"
}
stime<-mod$estimates
durata<-mod$time
#maxiter<-mod$maxiter
loglik<-as.numeric(mod$loglik)
niter<-mod$niter
varmat<-mod$varmat
BIC<-as.numeric(mod$BIC)
time<-mod$durata
results<-list('estimates'=stime,'ordinal'=ordinal,'time'=durata,
'loglik'=loglik,'niter'=niter,'varmat'=varmat,
'BIC'=BIC,'ellipsis'=ellipsis.arg,'family'=modelname,
'formula'=Formula,'call'=call)
attr(results,"hidden")<-c("estimates","ordinal","loglik","varmat","BIC","ellipsis","family")
class(results)<-"GEM"
# print.GEM<-function(x,...){
# hid<-attr(x,"hidden")
# print(x[!names(x)%in%hid])
# }
#
if (family == "cub"){
class(results)<-append("GEM","CUB")
#class(results)<-c("CUB","GEM")
}
if (family =="cube"){
class(results)<-append("GEM","CUBE")
}
if (family == "ihg"){
class(results)<-append("GEM","IHG")
}
if (family == "cush"){
class(results)<-append("GEM","CUSH")
}
invisible(results)
#return(results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.