Nothing
#' @title Main function for CUBE models with covariates only for feeling
#' @description Estimate and validate a CUBE model for ordinal data, with covariates only for explaining the
#' feeling component.
#' @aliases cubecsi
#' @usage cubecsi(m, ordinal, W, starting, maxiter, toler)
#' @param m Number of ordinal categories
#' @param ordinal Vector of ordinal responses
#' @param W Matrix of selected covariates for explaining the feeling component
#' @param starting Vector of initial parameters estimates to start the optimization algorithm, with length equal to
#' NCOL(W) + 3 to account for an intercept term for the feeling component (first entry)
#' @param maxiter Maximum number of iterations allowed for running the optimization algorithm
#' @param toler Fixed error tolerance for final estimates
#' @return An object of the class "CUBE". For cubecsi, $niter will return a NULL value since the optimization procedure
#' is not iterative but based on "optim" (method = "L-BFGS-B", option hessian=TRUE). \cr $varmat will return the inverse
#' of the numerically computed Hessian when it is positive definite, otherwise the procedure will return a matrix of NA
#' entries.
#' @import stats
#' @seealso \code{\link{loglikcubecsi}}, \code{\link{inibestcubecsi}}, \code{\link{CUBE}}
#' @keywords internal #models
cubecsi<-function(m,ordinal,W,starting,maxiter,toler){
tt0<-proc.time()
n<-length(ordinal)
W<-as.matrix(W)
if (ncol(W)==1){
W<-as.numeric(W)
}
q<-length(starting)-3
pai<-starting[1]; gama<-starting[2:(q+2)]; phi<-starting[q+3];
#(0)# log-lik
loglikzero<-loglikcubecsi(m,ordinal,W,pai,gama,phi)
#################################################################
param<-c(pai,gama,phi)
##################################
### maximize w.r.t. gama and phi #########
optimparam<-optim(param,effecubecsi,ordinal=ordinal,W=W,m=m,method="L-BFGS-B",lower=c(0.01,rep(-Inf,q+1),0.01), upper=c(0.99,rep(Inf,q+1),0.3),gr=NULL,hessian=TRUE)
#################################################################
# 7.# Computation of updated estimates and log-likelihood
#################################################################
paramest<-optimparam$par
pai<-paramest[1]
gama<-paramest[2:(q+2)] #updated gama estimates
phi<-paramest[q+3] #updated phi estimates
hessian<-optimparam$hessian
### updated log-likelihood
loglik<-loglikcubecsi(m,ordinal,W,pai,gama,phi)
vettestim<-c(pai,gama,phi)
nparam<-length(vettestim)
####################################################################
AICCUBEcsi<- -2*loglik+2*nparam
BICCUBEcsi<- -2*loglik+log(n)*nparam
###########################################################################
# Compute asymptotic standard errors of ML estimates via (numerical)Hessian
###########################################################################
if (det(hessian)<=0){
warning("Variance-Covariance matrix is not positive definite")
varmat<-ddd<-cormat<-matrix(NA,nrow=nparam,ncol=nparam)
errstd<-wald<-pval<-rep(NA,nparam)
ICOMP<-trvarmat<-NA
} else {
varmat<-solve(hessian)
errstd<-sqrt(diag(varmat))
ddd<-diag(sqrt(1/diag(varmat)))
wald<-vettestim/errstd
pval<-2*(1-pnorm(abs(wald)))
cormat<-(ddd%*%varmat)%*%ddd
trvarmat<-sum(diag(varmat))
ICOMP<- -2*loglik + nparam*log(trvarmat/nparam) - log(det(varmat))
errstd<-errstd
wald<-wald
pval<-pval
}
# nomi<-c("pai ",paste("gamma",0:(length(gama)-1),sep="_"),"phi ")
stime<-vettestim
####################################################################
### Print CUBEcsi results of ML estimation
####################################################################
# rownames(cormat)<-nomi; colnames(cormat)<-nomi;
durata<-proc.time()-tt0;durata<-durata[1];
results<-list('estimates'=stime, 'loglik'=loglik, 'varmat'=varmat,
'BIC'= BICCUBEcsi,'time'=durata,'niter'=1)
}
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.