############################################################################# ### >>> CUB models INFERENCE <<< Version 3.0 = September 2013 ### ############################################################################# #--------------------------------------------------------------------------- #(C) Copyright by Maria Iannario and Domenico Piccolo: (c)2009-2013 # maria.iannario@unina.it; domenico.piccolo@unina.it #*************************************************************************** #** All Rights Reserved. # (C) #** Version 1.0................February 2005 #** Version 2.0................September 2009 #** Version 3.0................September 2013 #** =================================================================================== #** This software derives from original GAUSS code, implemented by D.Piccolo(2001-2007) #** Shelter code has been implemented by M. Iannario...........March 2008-July 2009 #** last version (3.0) has been implemented by M. Iannario and D. Piccolo...Sept.2013 #** =================================================================================== #*************************************************************************** # !!!!!!!!!!!!!! Last modification on 31 August 2013 !!!!!!!!!!!!! #*************************************************************************** #** This entire File Header must accompany all files using #** any portion, in whole or in part, of this Source Code. #** If you wish to distribute any portion of the proprietary Source #** Code, you must first obtain written permission from the Author. #** #** If you use this program or part of it for applied or research projects, #** please inform the Author via E-mail to allow updating references file. #** When using this software, please correctly acknowledge it. #** #** This program is freely released without any responsibility for Authors. ############################################################################ #** Purpose: Estimate and test several models for ordinal data belonging # to the family of CUB models, also with shelter effect ############################################################################ #** ================================================================================ # The main reference for a correct usage of this program is: # > M.Iannario and D. Piccolo (2009), "A program in R for CUB models inference", # (version 2.0) # > M.Iannario and D. Piccolo (2013), "A Short Guide to CUB 3.0 Program". # #******************************************************************************* #**.... This program is currently maintained by: ..............................# #**.... M.Iannario ==> maria.iannario@unina.it; ...............................# #**.... D.Piccolo ==> domenico.piccolo@unina.it ..............................# #******************************************************************************* # We gratefully acknowledge people who positively improved previous versions #******************************************************************************* #** Notice that only the two global main functions (CUB and CUBE) are in upper case. #** ............................................................................... #** In general, ordinal responses are ratings in {1,2,...,m}, m>3. #** If examined as marginal values, univariate marginal ranking analyses are also possible #** (in this case, please be careful with interpretation). #** #** We denote by ordinal = (r_1,r_2,...,r_n)' the ordered (sampled) responses. #** #** If reverse rating/ranking is necessary, put: ordinal = m - ordinal + 1 #** .................................................................. #** If only frequencies=c(n_1,n_2,...,n_m) are available, then use: #** > ordinal=rep(1:m,frequencies); CUB(ordinal); #*********************************************************************************** #** Main notation adopted (for more details, see: Iannario and Piccolo, 2013): #** .................................................................. #** ......Y(n,p)=explanatory matrix variables for pai parameter......... #** ......W(n,q)=explanatory matrix variables for csi parameter......... #** ......X(n,s)=explanatory matrix variables for phi parameter......... #** .................................................................... #** .... Y, W, X do not require a vector of 1's in the first column !!! #** .................................................................... #** ....shelter=c (category, in {1,2,...,m} for a possible shelter effect) #** .................................................................... #** ... If Y or W are single dichotomous (=0,1) covariates, the software #** .... automatically produces more information and plots.............. #** .................................................................... #** ....Additional functions are available for advanced plots and simulation.... #** ================================================================================= #** > source("CUB.R") activates the program in R #** ================================================================================= #** Globals: m>3 scalar, number of categories, items, degrees, levels, etc. #** maxiter scalar, maximum number of E-M iterations. #** toler scalar, tolerance for relative change in parameters. ##################################################################################### # !!! m is a global variable that must be always declared !!! #### ##################################################################################### # *** Defaults: toler = 10e-06; maxiter = 500; ***** ##################################################################################### #** Main USAGE (without options):.................................................. #** .................................................................... #** In any case:==> > m=number-of-categorie; source("CUB.R"); #** CUB(0,0) ==> > CUB(ordinal); #** CUB(p,0) ==> > CUB(ordinal,Y=paicov); #** CUB(0,q) ==> > CUB(ordinal,W=csicov); #** CUB(p,q) ==> > CUB(ordinal,Y=paicov,W=csicov); #** CUB+shelter ==> > CUB(ordinal,shelter=c); #** CUBE(0,0,0) ==> > CUBE(ordinal); #** CUBE with covariates> CUBE(ordinal,Y=paicov,W=csicov,X=phicov); #** .................................................................... #============================================================================== #============================================================================== # *********************************************************************** # *** ................ SRUCTURE OF THE PROGRAM ...................... *** # *********************************************************************** #** I. General functions and statistical indices.................... #** II. Probability distributions.................................... #** III. Log-Likelihood functions..................................... #** IV. Var-covariance computations.................................. #** V. Initial values of estimates.................................. #** VI. Plotting facilities.......................................... #** VII Simulation routines.......................................... #** VIII. Main CUBs functions.......................................... #** IX Main call for CUB and CUBE models............................ #============================================================================== #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part I: General functions and statistical indices #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### DIRPROD ==> Direct product of a matrix with a vector dirprod<-function(Amat,xvett){ ra=nrow(Amat);ca=ncol(Amat); dprod=matrix(0,ncol=ca,nrow=ra); for(i in 1:ra){dprod[i,]=Amat[i,1:ca]*xvett[i]}; return(dprod); } ### LOGIS ==> Logistic function (...gives a vector whose lenght is rows(Y)..) logis<-function(Y,param){ YY=cbind(1,Y) # add 1's first column to matrix Y 1/(1+exp(-YY%*%param)) } ### KKK ==> Sequence of combinatorial coefficients kkk<-function(m,ordinal){choose(m-1,ordinal-1)} ### BETAR ==> Beta-Binomial probability distribution betar<-function(m,csi,phi){ betar=rep(NA,m); km=0:(m-2); betar[1]=prod(1-(1-csi)/(1+phi*km)); for(r in 1:(m-1)){ betar[r+1]=betar[r]*((m-r)/r)*((1-csi+phi*(r-1))/(csi+phi*(m-r-1))); } return(betar) } ### BETABINOMIAL... (vector) Beta-Binomial of ordinal, given csivett, phivett betabinomial<-function(m,csivett,phivett,ordinal){ n=length(ordinal); betabin=c(); for(i in 1:n){ bebeta=betar(m,csivett[i],phivett[i]) betabin[i]=bebeta[ordinal[i]] ### Attention !!! } return(betabin) } ### GINI ==> Normalized Gini (heterogeneity) index gini<-function(prob){(1-sum(prob^2))*m/(m-1)} ### LAAKSO ==> Normalized Laakso & Taagepera (heterogeneity) indexes laakso<-function(prob){1/(m/gini(prob)-m+1)} ### EXPCUB00 ==> Expectation of cub00 expcub00<-function(m,pai,csi){(m-1)*pai*(0.5-csi)+(m+1)/2} ### EXPCUBE==> Expectation of a CUBE model expcube<-function(m,pai,csi,phi){(m-1)*pai*(0.5-csi)+(m+1)/2} ### VARCUB00 ==> Varianceof cub00 varcub00<-function(m,pai,csi){ (m-1)*(pai*csi*(1-csi)+(1-pai)*(m+1)/12+pai*(1-pai)*(m-1)*((1-2*csi)^2)/4) } ### DISSOM ==> Normalized dissimilarity index dissom<-function(proba,probb){0.5*sum(abs(proba-probb))} ### AUP ==> Area Under Profile based on modal value (mode) for aggregated data ### observ=absolute frequencies; prob=theoretical probability aup<-function(m,observ,prob){ value=which.max(prob); ordinal=rep(1:m,observ); hitprof=c(); hitprof[1]=mean(abs(ordinal-value)==0); for(j in 1:m){ hitprof[j+1]=mean(abs(ordinal-value)<=j) } ### print(hitprof); (1+hitprof[1]+2*sum(hitprof[2:m]))/m-1 } #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part II: Probability distributions #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== #******************************************************** #** CUB ** probcub00=pai*probbit(m,csi)+(1-pai)/m *** #** CUBE ** probcube=pai*betar(csi,phi)+(1-pai)/m *** #******************************************************** ### PROBIT ==> Probability distribution of a shifted Binomial probbit<-function(m,csi){dbinom(0:(m-1),m-1,1-csi)} ### PROBCUB00 ==> Probability distribution of a CUB(0,0) model probcub00<-function(m,pai,csi){pai*(probbit(m,csi)-1/m)+1/m} ### BITCSI ==> Shifted Binomial with ordinal sequences bitcsi<-function(m,ordinal,csi){ base=log(1-csi)-log(csi); const=exp(m*log(csi)-log(1-csi)); const*kkk(m,ordinal)*exp(base*ordinal); } ### PROBCUBp0 ==> Probability distribution of a CUB(p,0) model probcubp0<-function(m,ordinal,Y,bet,csi){ logis(Y,bet)*(bitcsi(m,ordinal,csi)-1/m)+1/m; } ### BITGAMA ==> Shifted Binomial with ordinal sequences bitgama<-function(ordinal,W,gama){ WW=cbind(1,W); # add 1's first column to W ci=exp(-WW%*%gama); kkk(m,ordinal)*exp((ordinal-1)*log(ci)-(m-1)*log(1+ci)); } ### PROBCUB0q ==> Probability distribution of a CUB(0,q) model probcub0q<-function(m,ordinal,W,pai,gama){ pai*(bitgama(ordinal,W,gama)-1/m)+1/m; } ### PROBCUBpq ==> Probability distribution of a CUB(p,q) model probcubpq<-function(m,ordinal,Y,W,bet,gama){ logis(Y,bet)*(bitgama(ordinal,W,gama)-1/m)+1/m; } ### PROBCUBSHE1 ==> Probability distribution (1) of a CUB model + shelter effect ### see: Iannario M. (2012), SMA, 22, pp.5-7. probcubshe1<-function(m,pai1,pai2,csi,ccc) {pai1*probbit(m,csi)+pai2*(1/m)+(1-pai1-pai2)*ifelse(seq(1,m)==ccc,1,0)} ### PROBCUBSHE2 ==> Probability distribution (2) of a CUB model + shelter effect ### see: Iannario M. (2012), SMA, 22, pp.5-7. probcubshe2<-function(m,pai,csi,delta,shelter) {delta*ifelse(seq(1,m)==shelter,1,0)+(1-delta)*(pai*probbit(m,csi)+(1-pai)*(1/m))} ### PROBCUBE ==> Probability distribution of a CUBE model probcube<-function(m,pai,csi,phi){ pai*(betar(m,csi,phi)-1/m)+1/m } #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part III: Log-Likelihood functions #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### LOGLIKCUB00 ==> Log-likelihood function of a CUB(0,0) model loglikcub00<-function(freq,pai,csi){t(freq)%*%log(probcub00(m,pai,csi))} ### ELLECUB ==> Log-likelihood function of a CUB(0,0) model ### when ordinal data are available as r_i, i=1,2,...,n ellecub<-function(assepai,assecsi,ordinal){ prob=assepai*(dbinom(0:(m-1),m-1,1-assecsi)-1/m)+1/m; pconi=prob[ordinal]; sum(log(pconi)) } ### LOGLIKCUBp0 ==> Log-likelihood function of a CUB(p,0) model loglikcubp0<-function(ordinal,Y,bbet,ccsi){ prob=probbit(m,ccsi); probn=prob[ordinal]; YY=cbind(1,Y); # add 1's first column to Y eta=1/(1+exp(-YY%*%bbet)); sum(log(eta*(probn-1/m)+1/m)); } ### LOGLIKCUB0q ==> Log-likelihood function of a CUB(0,q) model loglikcub0q<-function(ordinal,W,pai,gama){ probn=probcub0q(m,ordinal,W,pai,gama); sum(log(probn)); } ### LOGLIKCUBpq ==> Log-likelihood function of a CUB(p,q) model loglikcubpq<-function(ordinal,Y,W,bet,gama){ probn=probcubpq(m,ordinal,Y,W,bet,gama); sum(log(probn)); } ### LOGLIKCUBSHE ==> Log-likelihood function of a CUB model + shelter effect loglikcubshe<-function(freq,pai1,pai2,csi,ccc) {t(freq)%*%log(probcubshe1(m,pai1,pai2,csi,ccc))} ### LOGLIKCUBE ==> Log-likelihood function for frequencies: (freq_i, i=1,2,...,n) loglikcube<-function(m,pai,csi,phi,freq){t(freq)%*%log(probcube(m,pai,csi,phi))} ### LOGLIKECUBEN ==> Log-likelihood function for ordinal data: (r_i, i=1,2,...,n) loglikcuben<-function(m,assepai,assecsi,assephi,ordinal){ prob=probcube(m,assepai,assecsi,assephi); pconi=prob[ordinal]; sum(log(pconi)) } #### EFFECUBE ==> Q(csi,phi) ### dati=cbind(tauno,freq) ############# EFFECUBE<-function(paravec,dati){ tauno=dati[,1]; freq=dati[,2]; -sum(freq*tauno*log(betar(m,paravec[1],paravec[2]))) } ### LOGLIKECUBECOV ==> (scalar) Log-likelihood function of CUBE models ### for ordinal with (Y,W,Z) loglikcubecov<-function(ordinal,Y,W,Z,bet,gama,alpha){ paivett=logis(Y,bet); csivett=logis(W,gama); phivett=logis(Z,alpha); probi=paivett*(betabinomial(m,csivett,phivett,ordinal)-1/m)+1/m; sum(log(probi)); } #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part IV: Var-covariance computations #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### VARCOVCUB00 ==> Variance-covariance matrix for estimation of a CUB(0,0) model varcovcub00<-function(ordinal,pai,csi){ vvi=(m-ordinal)/csi-(ordinal-1)/(1-csi); ui=(m-ordinal)/(csi^2)+(ordinal-1)/((1-csi)^2); pri=probcub00(m,pai,csi); qi=1/(m*pri[ordinal]); qistar=1-(1-pai)*qi; qitilde=qistar*(1-qistar); i11=sum((1-qi)^2)/(pai^2); i12=-sum(vvi*qi*qistar)/pai; i22=sum(qistar*ui-(vvi^2)*qitilde); ####################################### Information matrix matinf=matrix(c(i11,i12,i12,i22),nrow=2,byrow=T); ####################################### Variance-covariance matrix notpd=0; if(det(matinf)<=0){ notpd=1; cat("=======================================================================","\n") cat("Variance-covariance matrix NOT positive definite","\n") cat("=======================================================================","\n") } else { notpd=0; varmat=solve(matinf); } } ### VARCOVCUBp0 ==> Variance-covariance matrix for estimation of a CUB(p,0) model varcovcubp0<-function(ordinal,Y,bet,csi){ vvi=(m-ordinal)/csi-(ordinal-1)/(1-csi); ui=(m-ordinal)/(csi^2)+(ordinal-1)/((1-csi)^2); qi=1/(m*probcubp0(m,ordinal,Y,bet,csi)); ei=logis(Y,bet); qistar=1-(1-ei)*qi; eitilde=ei*(1-ei); qitilde=qistar*(1-qistar); ff=eitilde-qitilde; g10=vvi*qitilde; YY=cbind(1,Y); ### add 1's first column to Y i11=t(YY)%*%(dirprod(YY,ff)); ### i11=t(YY)%*%(YY*ff); i12=-t(YY)%*%(g10); i22=sum(ui*qistar-(vvi^2)*qitilde); # Information matrix matinf=rbind(cbind(i11,i12),cbind(t(i12),i22)); # Var-covar matrix varmat=solve(matinf); return(varmat); } ### VARCOVCUB0q ==> Variance-covariance matrix for estimation of a CUB(0,q) model varcovcub0q<-function(ordinal,W,pai,gama){ qi=1/(m*probcub0q(m,ordinal,W,pai,gama)); qistar=1-(1-pai)*qi; qitilde=qistar*(1-qistar); fi=logis(W,gama); fitilde=fi*(1-fi); ai=(ordinal-1)-(m-1)*(1-fi); g01=(ai*qi*qistar)/pai; hh=(m-1)*qistar*fitilde-(ai^2)*qitilde; WW=cbind(1,W); # add 1's first column to W i11=sum((1-qi)^2)/(pai^2); i12=t(g01)%*%WW; i22=t(WW)%*%(dirprod(WW,hh)); ### i22=t(WW)%*%(WW*hh); matinf=rbind(cbind(i11,i12),cbind(t(i12),i22)); ### Information matrix varmat=solve(matinf); ### Var-covar matrix return(varmat); } ### VARCOVCUBpq ==> Variance-covariance matrix for estimation of a CUB(p,q) model varcovcubpq<-function(ordinal,Y,W,bet,gama){ qi=1/(m*probcubpq(m,ordinal,Y,W,bet,gama)); ei=logis(Y,bet); eitilde=ei*(1-ei); qistar=1-(1-ei)*qi; qitilde=qistar*(1-qistar); fi=logis(W,gama); fitilde=fi*(1-fi); ai=(ordinal-1)-(m-1)*(1-fi); ff=eitilde-qitilde; gg=ai*qitilde; hh=(m-1)*qistar*fitilde-(ai^2)*qitilde; YY=cbind(1,Y); # add 1's first column to Y WW=cbind(1,W); # add 1's first column to W i11=t(YY)%*%(dirprod(YY,ff)); ### i11=t(YY)%*%(YY*ff); i12=t(YY)%*%(dirprod(WW,gg)); ### i12=t(YY)%*%(WW*gg); i22=t(WW)%*%(dirprod(WW,hh)); ### i22=t(WW)%*%(WW*hh); matinf=rbind(cbind(i11,i12),cbind(t(i12),i22)); # Information matrix varmat=solve(matinf); # Var-covar matrix return(varmat); } ### VARCOVCUBSHE ==> Variance-covariance matrix for estimation of a CUB model + shelter effect varcovcubshe<-function(pai1,pai2,csi,ccc,n){ pr=probcubshe1(m,pai1,pai2,csi,ccc); dd=rep(0,m);dd[ccc]=1; bb=probbit(m,csi); ######################## aaa=bb-dd; bbb=(1/m)-dd; c4=pai1*bb*(m-(1:m)-csi*(m-1))/(csi*(1-csi)); atilde=aaa/pr; btilde=bbb/pr; ctilde=c4/pr; d11=sum(aaa*atilde); d22=sum(bbb*btilde); dxx=sum(c4*ctilde); d12=sum(bbb*atilde); d1x=sum(c4*atilde); d2x=sum(c4*btilde); ### Information matrix matinf=matrix(c(d11,d12,d1x,d12,d22,d2x,d1x,d2x,dxx),nrow=3,byrow=T); ### Var-covar matrix if(det(matinf)<=0){ notpd=1; cat("=======================================================================","\n") cat("Variance-covariance matrix NOT positive definite","\n") cat("=======================================================================","\n") assign('notpd',notpd,pos=1); } else { notpd=0; } ### Attention!!! 1/n varmat=solve(matinf)/n; return(varmat) } ################################################################ ### EXPECTED Information matrix (optional computation) ### VARCOVCUBEEXP ==> var-covar from EXPECTED information matrix ### for a CUBE(pai,csi,phi) on n ratings ################################################################ varcovcubeexp<-function(m,pai,csi,phi,n){ pr=probcube(m,pai,csi,phi); sum1=sum2=sum3=sum4=c(); for(jr in 1:m){ seq1=1/((1-csi)+phi*((1:jr)-1)); seq2=1/((csi)+phi*((1:(m-jr+1))-1)); seq3=((1:jr)-1)/((1-csi)+phi*((1:jr)-1)); seq4=((1:(m-jr+1))-1)/((csi)+phi*((1:(m-jr+1))-1)); sum1[jr]=sum(seq1)-seq1[jr]; sum2[jr]=sum(seq2)-seq2[m-jr+1]; sum3[jr]=sum(seq3)-seq3[jr]; sum4[jr]=sum(seq4)-seq4[m-jr+1]; } sum5=sum(((1:m)-1)/(1+phi*(1:m))); cr=pr-(1-pai)/m; derpai=(pr-1/m)/pai; dercsi=cr*(sum2-sum1); derphi=cr*(sum3+sum4-sum5); infpaipai=sum((derpai^2)/pr); infpaicsi=sum((derpai*dercsi)/pr); infpaiphi=sum((derpai*derphi)/pr); infcsicsi=sum((dercsi^2)/pr); infcsiphi=sum((dercsi*derphi)/pr); infphiphi=sum((derphi^2)/pr); ### Information matrix inform=matrix(NA,nrow=3,ncol=3); inform[1,1]=infpaipai; inform[1,2]=infpaicsi; inform[1,3]=infpaiphi; inform[2,1]=infpaicsi; inform[2,2]=infcsicsi; inform[2,3]=infcsiphi; inform[3,1]=infpaiphi; inform[3,2]=infcsiphi; inform[3,3]=infphiphi; ### Var-covar matrix varmat=(solve(inform))/n; return(varmat) } ############################################################# ### OBSERVED Information Matrix (default computation) ############################################################# ### VARCOVCUBEOBS ==> var-covar from OBSERVED information matrix ### for a CUBE(pai,csi,phi) on n ratings ############################################################# varcovcubeobs<-function(m,pai,csi,phi,freq){ pr=probcube(m,pai,csi,phi); sum1=sum2=sum3=sum4=c(); ### sums computation d1=d2=h1=h2=h3=h4=c(); ### derivatives computation ### sum1; sum2; sum3; sum4; sum5; as in Iannario (2013), "Comm. in Stat. Theory & Methods" for(jr in 1:m){ seq1=1/((1-csi)+phi*((1:jr)-1)); ### a(k) seq2=1/((csi)+phi*((1:(m-jr+1))-1)); ### b(k) seq3=((1:jr)-1)/((1-csi)+phi*((1:jr)-1)); seq4=((1:(m-jr+1))-1)/((csi)+phi*((1:(m-jr+1))-1)); dseq1=seq1^2; dseq2=seq2^2; hseq1=dseq1*((1:jr)-1); hseq2=dseq2*((1:(m-jr+1))-1); hseq3=dseq1*((1:jr)-1)^2; hseq4=dseq2*((1:(m-jr+1))-1)^2; ############# sum1[jr]=sum(seq1)-seq1[jr]; sum2[jr]=sum(seq2)-seq2[m-jr+1]; sum3[jr]=sum(seq3)-seq3[jr]; sum4[jr]=sum(seq4)-seq4[m-jr+1]; d1[jr]= sum(dseq1)-dseq1[jr]; d2[jr]=-(sum(dseq2)-dseq2[m-jr+1]); h1[jr]=-(sum(hseq1)-hseq1[jr]); h2[jr]=-(sum(hseq2)-hseq2[m-jr+1]); h3[jr]=-(sum(hseq3)-hseq3[jr]); h4[jr]=-(sum(hseq4)-hseq4[m-jr+1]); } seq5=(0:(m-2))/(1+phi*(0:(m-2))); ### c(k) correction ?! sum5=sum(seq5); ### correction ?! h5=-sum(seq5^2); ### correction ?! ### Symbols as in Iannario (2013), "Comm. in Stat.", ibidem (DP notes) uuur=1-1/(m*pr) ubar=uuur+pai*(1-uuur) vbar=ubar-1 aaar=sum2-sum1 bbbr=sum3+sum4-sum5 cccr=h3+h4-h5 dddr=h2-h1 eeer=d2-d1 ###### dummy product prodo=freq*ubar ###### infpaipai=sum(freq*uuur^2)/pai^2; infpaicsi=sum(prodo*(uuur-1)*aaar)/pai; infpaiphi=sum(prodo*(uuur-1)*bbbr)/pai; infcsicsi=sum(prodo*(vbar*aaar^2-eeer)); infcsiphi=sum(prodo*(vbar*aaar*bbbr-dddr)); infphiphi=sum(prodo*(vbar*bbbr^2-cccr)); ### Information matrix inform=matrix(NA,nrow=3,ncol=3); inform[1,1]=infpaipai; inform[1,2]=infpaicsi; inform[1,3]=infpaiphi; inform[2,1]=infpaicsi; inform[2,2]=infcsicsi; inform[2,3]=infcsiphi; inform[3,1]=infpaiphi; inform[3,2]=infcsiphi; inform[3,3]=infphiphi; ### Var-covar matrix varmat=solve(inform); return(varmat) } ### ELEMENTO ==> Produce a matrix(m,n) whose element is: ### matrix[k,i]=(e*(k-1)^d)/(a+b*csi_i+phi_i*(k-1))^c elemento<-function(m,vettcsi,vettphi,a,b,c,d,e){ elemat=matrix(NA,nrow=m,ncol=length(vettcsi)) for(k in 1:m){ elemat[k,]=e*((k-1)^d)/((a+b*vettcsi+vettphi*(k-1))^c) } return(elemat) } ### VARCOVCUBECOV ==> Compute CUBE var-cov matrix of estimates varcovcubecov<-function(m,ordinal,Y,W,Z,estbet,estgama,estalpha){ n=length(ordinal); paivett=logis(Y,estbet); csivett=logis(W,estgama); phivett=logis(Z,estalpha); ### Probability probi=paivett*(betabinomial(m,csivett,phivett,ordinal)-1/m)+1/m; uui=1-1/(m*probi) ubari=uui+paivett*(1-uui) ### Matrix computations mats1=elemento(m,csivett,phivett,1,-1,1,0,1) mats2=elemento(m,csivett,phivett,0, 1,1,0,1) mats3=elemento(m,csivett,phivett,1,-1,1,1,1) mats4=elemento(m,csivett,phivett,0, 1,1,1,1) mats5=elemento(m,csivett,phivett,1, 0,1,1,1) ### per vettori D_i matd1=elemento(m,csivett,phivett,1,-1,2,0, 1) matd2=elemento(m,csivett,phivett,0, 1,2,0,-1) matd3=elemento(m,csivett,phivett,1,-1,2,1, 1) matd4=elemento(m,csivett,phivett,0, 1,2,1,-1) ### per vettori H_i math3=elemento(m,csivett,phivett,1,-1,2,2,-1) math4=elemento(m,csivett,phivett,0, 1,2,2,-1) math5=elemento(m,csivett,phivett,1, 0,2,2,-1) ### S1=S2=S3=S4=D1=D2=D3=D4=H3=H4=c(); for(i in 1:length(ordinal)){ S1[i]=sum(mats1[1:ordinal[i],i])-mats1[ordinal[i],i]; S2[i]=sum(mats2[1:(m-ordinal[i]+1),i])-mats2[(m-ordinal[i]+1),i]; S3[i]=sum(mats3[1:ordinal[i],i])-mats3[ordinal[i],i]; S4[i]=sum(mats4[1:(m-ordinal[i]+1),i])-mats4[(m-ordinal[i]+1),i]; D1[i]=sum(matd1[1:ordinal[i],i])-matd1[ordinal[i],i]; D2[i]=sum(matd2[1:(m-ordinal[i]+1),i])-matd2[(m-ordinal[i]+1),i]; D3[i]=sum(matd3[1:ordinal[i],i])-matd3[ordinal[i],i]; D4[i]=sum(matd4[1:(m-ordinal[i]+1),i])-matd4[(m-ordinal[i]+1),i]; H3[i]=sum(math3[1:ordinal[i],i])-math3[ordinal[i],i]; H4[i]=sum(math4[1:(m-ordinal[i]+1),i])-math4[(m-ordinal[i]+1),i]; } ###### Attention !!! S5=colSums(mats5[1:(m-1),]); H5=colSums(math5[1:(m-1),]); ########## CC=S2-S1; EE=S3+S4-S5; DD=D2-D1; FF=D3+D4; GG=H3+H4-H5; ### Computing vectors v and u vibe=uui*(1-paivett); viga=ubari*csivett*(1-csivett)*CC; vial=ubari*phivett*(1-phivett)*EE; ubebe=uui*(1-paivett)*(1-2*paivett); ugabe=ubari*csivett*(1-csivett)*(1-paivett)*CC; ualbe=ubari*phivett*(1-phivett)*(1-paivett)*EE; ugaga=ubari*csivett*(1-csivett)* ((1-2*csivett)*CC+csivett*(1-csivett)*(CC^2+DD)); ualga=ubari*phivett*(1-phivett)*csivett*(1-csivett)*(FF+CC*EE); ualal=ubari*phivett*(1-phivett)* ((1-2*phivett)*EE+phivett*(1-phivett)*(EE^2+GG)); ### Computing vectors g gbebe=dirprod(vibe,vibe)-ubebe; ggabe=dirprod(viga,vibe)-ugabe; galbe=dirprod(vial,vibe)-ualbe; ggaga=dirprod(viga,viga)-ugaga; galga=dirprod(vial,viga)-ualga; galal=dirprod(vial,vial)-ualal; ### Expanding matrices Y, W, Z YY=cbind(1,Y); WW=cbind(1,W); ZZ=cbind(1,Z); ### Elements of the Information matrix infbebe=t(YY)%*%(dirprod(YY,gbebe)); infgabe=t(WW)%*%(dirprod(YY,ggabe)); infalbe=t(ZZ)%*%(dirprod(YY,galbe)); infgaga=t(WW)%*%(dirprod(WW,ggaga)); infalga=t(ZZ)%*%(dirprod(WW,galga)); infalal=t(ZZ)%*%(dirprod(ZZ,galal)); infbega=t(infgabe); infbeal=t(infalbe); infgaal=t(infalga); ### Assembling Information matrix... matinf=rbind( cbind(infbebe,infbega,infbeal), cbind(infgabe,infgaga,infgaal), cbind(infalbe,infalga,infalal)); ### Variance-covariance matrix varmat=solve(matinf); return(varmat) } #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part V: Initial values of estimates #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### INIBEST ==> Initial values of EM algorithm for a CUB(0,0) model inibest<-function(freq){ m=length(freq); freq=freq/sum(freq); csi=1+(0.5-which.max(freq))/m; ppp=probbit(m,csi); pai=sqrt((sum(freq^2)-1/m)/(sum(ppp^2)-1/m)); pai=min(pai,0.99); return(c(pai,csi)); } ### INIGRID ==> Grid-based initial values of EM algorithm for a CUB(0,0) model inigrid<-function(freq){ listap=list(c(0.1,0.1),c(0.1,0.25),c(0.1,0.5),c(0.1,0.65),c(0.1,0.85), c(0.25,0.1),c(0.25,0.25),c(0.25,0.5),c(0.25,0.65),c(0.25,0.85), c(0.5,0.1),c(0.5,0.25),c(0.5,0.5),c(0.5,0.65),c(0.5,0.85), c(0.65,0.1),c(0.65,0.25),c(0.65,0.5),c(0.65,0.65),c(0.65,0.85), c(0.85,0.1),c(0.85,0.25),c(0.85,0.5),c(0.85,0.65),c(0.1,0.85)) quanti=length(listap) for(j in 1:quanti){ pai=listap[[j]][1]; csi=listap[[j]][2] loglik[j]=loglikcub00(freq,pai,csi) } indice=which.max(loglik) return(listap[[indice]]) } ### INIBESTGAMA ==> Initial values of EM algorithm for a CUB(0,q) model inibestgama<-function(ordinal,W){ WW=cbind(1,W); # add 1's first column to W ni=log((m-ordinal+0.5)/(ordinal-0.5)); gama=(solve(t(WW)%*%WW))%*%(t(WW)%*%ni); return(gama); } ### INIBESTCUBE ==> Initial values of EM algorithm for a CUBE model inibestcube<-function(ordinal){ out=cub00forsim(ordinal,maxiter=500,toler=1e-6) initial=c(out[1],out[2],0.1) return(initial); } #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part VI: Plotting facilities #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### CUBVISUAL ==> plot an estimated CUB model as a point in the parameter space cubvisual<-function(ordinal,labpoint){ cub00(ordinal,makeplot=FALSE) plot(1-pai,1-csi,main="CUB models parameter space",las=1,pch=19,cex=1.5, xlim=c(0,1),ylim=c(0,1),col="blue", xlab=expression(paste("Uncertainty ", (1-pi))), ylab=expression(paste("Feeling ", (1-xi)))); text(1-pai,1-csi,labels=labpoint,font=4,pos=1,offset=0.4,cex=0.8) } ### DICOPAI ==> compares and plots 2 CUB(p,0) distributions ### when the unique covariate for pai is dichotomous (=0,1) dicopai<-function(ordinal,Y,m,bet,csi){ pai0=1/(1+exp(-(bet[1]+bet[2]*0))); prob0=probcub00(m,pai0,csi); pai1=1/(1+exp(-(bet[1]+bet[2]*1))); prob1=probcub00(m,pai1,csi); maxpr=max(prob0,prob1); makeplot=TRUE ### Alternatively, makeplot=FALSE if(makeplot==TRUE){ plot(1:m,prob0,ylim=c(0.0,1.1*maxpr),cex.main=0.8,las=1, main="CUB distributions, given pai-covariate=0, 1",cex=2, xlab="",ylab="Prob(R|D=0) and Prob(R|D=1)",pch=1,lty=1,type="b"); lines(1:m,prob1,cex=2,pch=19,lty=2,type="b"); abline(h=0); } ### Expected moments given D=0,1 exp0=expcub00(m,pai0,csi); exp1=expcub00(m,pai1,csi); cubmode0=which.max(prob0); cubmode1=which.max(prob1); ### Sample averages and modal value, given D=0,1 ord0=ordinal[Y==0]; ord1=ordinal[Y==1]; n0=length(ord0); n1=length(ord1); aver0=mean(ord0); aver1=mean(ord1); obsmode0=which.max(tabulate(ordinal[Y==0])); obsmode1=which.max(tabulate(ordinal[Y==1])) cat("Samples and populations measures, given dichotomous covariate (D=0) and (D=1)","\n"); cat("-----------------------------------------------------------------------","\n"); cat("(D = 0)"," n = ", round(n0,digits=1), " pai_0=",round(pai0,digits=3)," csi=",round(csi,digits=3),"\n"); cat("............................","\n"); cat("Sample average =",round(aver0,digits=8)," Sample mode =",round(obsmode0,digits=1),"\n"); cat("CUB expectation =",round(exp0,digits=8), " CUB mode =",round(cubmode0,digits=1),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("(D = 1)"," n = ", round(n1,digits=1), " pai_1=",round(pai1,digits=3)," csi=",round(csi,digits=3),"\n"); cat("............................","\n"); cat("Sample average =",round(aver1,digits=8)," Sample mode =",round(obsmode1,digits=1),"\n"); cat("CUB expectation =",round(exp1,digits=8), " CUB mode =",round(cubmode1,digits=1),"\n"); cat("-----------------------------------------------------------------------","\n"); } ### DICOCSI ==> compares and plots 2 CUB(0,q) distributions ### when the unique covariate for csi is dichotomous (=0,1) dicocsi<-function(ordinal,W,m,pai,gama){ csi0=1/(1+exp(-(gama[1]+gama[2]*0))); prob0=probcub00(m,pai,csi0); csi1=1/(1+exp(-(gama[1]+gama[2]*1))); prob1=probcub00(m,pai,csi1); maxpr=max(prob0,prob1); makeplot=TRUE ### Alternatively, makeplot=FALSE if(makeplot==TRUE){ plot(1:m,prob0,ylim=c(0.0,1.1*maxpr),cex.main=0.8,las=1, main="CUB distributions, given csi-covariate=0, 1", cex=2,xlab="",ylab="Prob(R|D=0) and Prob(R|D=1)",pch=1,lty=1,type="b"); lines(1:m,prob1,cex=2,pch=19,lty=2,type="b"); abline(h=0); } ### Expected moments given D=0,1 exp0=expcub00(m,pai,csi0); exp1=expcub00(m,pai,csi1); cubmode0=which.max(prob0); cubmode1=which.max(prob1); ### Sample averages and modal value, given D=0,1 ord0=ordinal[W==0]; ord1=ordinal[W==1]; n0=length(ord0); n1=length(ord1); aver0=mean(ord0); aver1=mean(ord1); obsmode0=which.max(tabulate(ord0)); obsmode1=which.max(tabulate(ord1)) cat("Samples and populations measures, given dichotomous covariate (D=0) and (D=1)","\n"); cat("-----------------------------------------------------------------------","\n"); cat("(D = 0)"," n = ", round(n0,digits=1), " pai=",round(pai,digits=3)," csi_0=",round(csi0,digits=3),"\n"); cat("............................","\n"); cat("Sample average =",round(aver0,digits=8)," Sample mode =",round(obsmode0,digits=1),"\n"); cat("CUB expectation =",round(exp0,digits=8), " CUB mode =",round(cubmode0,digits=1),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("(D = 1)"," n = ", round(n1,digits=1), " pai=",round(pai,digits=3)," csi_1=",round(csi1,digits=3),"\n"); cat("............................","\n"); cat("Sample average =",round(aver1,digits=8)," Sample mode =",round(obsmode1,digits=1),"\n"); cat("CUB expectation =",round(exp1,digits=8), " CUB mode =",round(cubmode1,digits=1),"\n"); cat("-----------------------------------------------------------------------","\n"); } ### DICOPAICSI ==> compares and plots 2 CUB(p,q) distributions ### when the same covariate for pai and csi is dichotomous (=0,1) dicopaicsi<-function(ordinal,Y,W,m,bet,gama){ pai0=1/(1+exp(-(bet[1]+bet[2]*0))); pai1=1/(1+exp(-(bet[1]+bet[2]*1))); csi0=1/(1+exp(-(gama[1]+gama[2]*0))); csi1=1/(1+exp(-(gama[1]+gama[2]*1))); prob0=probcub00(m,pai0,csi0); prob1=probcub00(m,pai1,csi1); maxpr=max(prob0,prob1); makeplot=TRUE ### Alternatively, makeplot=FALSE if(makeplot==TRUE){ plot(1:m,prob0,ylim=c(0.0,1.1*maxpr),cex.main=0.8,las=1, main="CUB distributions, given pai and csi covariates=0, 1",cex=2, xlab="",ylab="Prob(R|D=0) and Prob(R|D=1)",pch=1,lty=1,type="b"); lines(1:m,prob1,cex=2,pch=19,lty=2,type="b"); abline(h=0); } ### Expected moments given D=0,1 exp0=expcub00(m,pai0,csi0); exp1=expcub00(m,pai1,csi1); cubmode0=which.max(prob0); cubmode1=which.max(prob1); ### Sample averages and modal value, given D=0,1 ord0=ordinal[Y==0]; ord1=ordinal[Y==1]; n0=length(ord0); n1=length(ord1); aver0=mean(ord0); aver1=mean(ord1); obsmode0=which.max(tabulate(ordinal[Y==0])); obsmode1=which.max(tabulate(ordinal[Y==1])) cat("Samples and populations measures, given dichotomous covariate (D=0) and (D=1)","\n"); cat("-----------------------------------------------------------------------","\n"); cat("(D = 0)"," n = ", round(n0,digits=1), " pai_0=",round(pai0,digits=3)," csi_0=",round(csi0,digits=3),"\n"); cat("............................","\n"); cat("Sample average =",round(aver0,digits=8)," Sample mode =",round(obsmode0,digits=1),"\n"); cat("CUB expectation =",round(exp0,digits=8), " CUB mode =",round(cubmode0,digits=1),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("(D = 1)"," n = ", round(n1,digits=1), " pai_1=",round(pai1,digits=3)," csi_1=",round(csi1,digits=3),"\n"); cat("............................","\n"); cat("Sample average =",round(aver1,digits=8)," Sample mode =",round(obsmode1,digits=1),"\n"); cat("CUB expectation =",round(exp1,digits=8), " CUB mode =",round(cubmode1,digits=1),"\n"); cat("-----------------------------------------------------------------------","\n"); } ### MULTICUB ==> Joint plotting of estimated CUB models in the parameter space multicub<-function(matord,m,etich=as.character(1:ncol(matord)), titolo="CUB models",colori="black", simboli=19, thickness=1.5,xwidth=c(0,1),ywidth=c(0,1)){ k=ncol(matord); vettpai=vettcsi=c() for(j in 1:k){ cub00(matord[,j],makeplot=FALSE) vettpai[j]=pai; vettcsi[j]=csi; } plot(1-vettpai,1-vettcsi,main=titolo,cex=1.2,cex.main=1, font.lab=4,cex.lab=1, pch=simboli, col=colori,lwd=thickness,las=1, xlim=xwidth,ylim=ywidth, xlab=expression(paste("Uncertainty ", (1-pi))), ylab=expression(paste("Feeling ", (1-xi)))); text(1-vettpai,1-vettcsi,labels=etich,pos=3,offset=0.5,font=4,cex=0.9) } ### MULTICUBE ==> Joint plotting of estimated CUBE models in the parameter space multicube<-function(matord,m,etich=as.character(1:ncol(matord)), titolo="CUBE models",colori="black",simboli=19, thickness=1.5,xwidth=c(0,1),ywidth=c(0,1)){ k=ncol(matord); vettpai=vettcsi=c() for(j in 1:k){ CUBE(matord[,j],makeplot=FALSE) vettpai[j]=pai; vettcsi[j]=csi; } plot(1-vettpai,1-vettcsi,main=titolo,cex=1.2,cex.main=1, font.lab=4,cex.lab=1, pch=simboli, col=colori,lwd=thickness,las=1, xlim=xwidth,ylim=ywidth, xlab=expression(paste("Uncertainty ", (1-pi))), ylab=expression(paste("Feeling ", (1-xi)))); text(1-vettpai,1-vettcsi,labels=etich,pos=2,offset=0.3,cex=0.7) } #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part VII: Simulation routines #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### SIMCUB ==> Simulation of n pseudo-random numbers generated ### by a CUB model without covariates simcub<-function(n,m,pai,csi){ dico=runif(n) Simulation of n pseudo-random numbers generated ### by a CUB model + shelter effect simcubshe<-function(n,m,pai,csi,delta,shelter){ dicopai=runif(n) Simulation of n pseudo-random numbers generated ### by a CUBE model simcube<-function(n,m,pai,csi,phi){ prob=probcube(m,pai,csi,phi); sample(1:m,n,prob,replace=T) } ### CUB00FORSIM ==> CUB Estimation routine, useful for simulation experiments cub00forsim<-function(ordinal,maxiter=500,toler=1e-6){ tt0=proc.time(); serie=1:m; freq=tabulate(ordinal,nbins=m); n=sum(freq); aver=mean(ordinal); varcamp=mean(ordinal^2)-aver^2; ############################################################## inipaicsi=inibest(freq); pai=inipaicsi[1]; csi=inipaicsi[2]; ### pai=runif(1);csi=runif(1); ==> Random initial values ################################################################## loglik=loglikcub00(freq,pai,csi); # ******************************************************************** # ************* E-M algorithm for CUB(0,0) *************************** # ******************************************************************** nniter=1; while(nniter<=maxiter){ likold=loglik; bb=probbit(m,csi); aa=(1-pai)/(m*pai*bb); tau=1/(1+aa); ft=freq*tau; averpo=(t(serie)%*%ft)/sum(ft); pai=(t(freq)%*%tau)/n; # updated pai estimate csi=(m-averpo)/(m-1); # updated csi estimate if(csi<0.001){ csi=0.01; nniter=maxiter-1 } # print(c(pai,csi)); loglik=loglikcub00(freq,pai,csi); liknew=loglik; testll=abs(liknew-likold); ###### print(testll); # OPTIONAL printing: print(cbind(nniter,testll,pai,csi)); if(testll<=toler) break else {loglik=liknew}; # OPTIONAL printing: print(loglik); nniter=nniter+1; } ###### if(csi>0.999) csi=0.99; ### to avoid division by 0 !!! if(csi<0.001) csi=0.01; ### to avoid division by 0 !!! if(pai<0.001) pai=0.01; ### to ensure identifiability !!! ################################################################ assign('loglik',loglik,pos=1); assign('pai',pai,pos=1); assign('csi',csi,pos=1); return(c(pai,csi)) } ###### >>>>> End of CUB00FORSIM <<<<< ####################### ### CUBEFORSIM ==> CUBE Estimation routine, useful for simulation experiments cubeforsim<-function(ordinal,starting=rep(0.1,3),maxiter=500,toler=1e-6){ freq=tabulate(ordinal,nbins=m); n=sum(freq); aver=mean(ordinal); varcamp=mean(ordinal^2)-aver^2; ######################################################## #* Defaults: toler = 0.000001; maxiter = 500; ***** ######################################################## #(0)# initial estimates pai=0.1; csi=0.1; phi=0.1; #(0)# log-lik loglik=loglikcube(m,pai,csi,phi,freq); # ******************************************************************** # ************* E-M algorithm for CUBE ************************* # ******************************************************************** nniter=1; while(nniter<=maxiter){ likold=loglik; #(1)# betar bb=betar(m,csi,phi); aa=(1-pai)/(m*pai*bb); #(2)# taunor tauno=1/(1+aa); #(3)# pai(k+1) pai=sum(freq*tauno)/n; # updated pai estimate paravecjj=c(csi,phi); #(4)# Q(k+1) dati=cbind(tauno,freq); ################ EFFECUBE is Q(csi,phi) ########################### #(5)# (csi(k+1),phi(k+1)) ################################## maximize w.r.t. paravec ######## optimestim=optim(paravecjj,EFFECUBE,dati=dati); # print(nlmaxg) ### prova simulazione ### ### optimestim=optim(paravecjj,EFFECUBE,dati=dati,method = "L-BFGS-B",lower=c(0.01,0.1),upper=c(0.99,0.3)) ################################################################ #(6)# theta(k+1) paravecjj=optimestim$par; # updated paravec estimates csi=paravecjj[1]; phi=paravecjj[2]; ########################################## # if(pai<0.001){pai=0.001; nniter=maxiter-1} # if(csi<0.001){csi=0.001; nniter=maxiter-1} # if(phi<0.001){phi=0.001; nniter=maxiter-1} # if(pai>0.999){pai=0.98} ### to avoid division by 0 !!! # if(csi>0.999){csi=0.98} ### to avoid division by 0 !!! ###################################### print(c(nniter,pai,csi,phi)) if(phi>1000){phi=1000; nniter=maxiter-1} ### to avoid not-convergence !!! #(7)# elle(theta(k+1)) liknew=loglikcube(m,pai,csi,phi,freq); #(8)# test testll=abs(liknew-likold); # OPTIONAL printing: print(testll); # OPTIONAL printing: print(cbind(nniter,testll,pai,csi,phi)); if(testll<=toler) break else {loglik=liknew}; # OPTIONAL printing: print(loglik); nniter=nniter+1; } loglik=liknew; ###### End of E-M algorithm for CUBE *********************************************** # ******************************************************** # Compute ML var-cov matrix and print result for CUBE # ******************************************************** varmat=try(varcovcubeobs(m,pai,csi,phi,freq)) ##################################################################### # Assignments as global variables: assign('name',value,pos=1); ##################################################################### assign('pai',pai,pos=1); assign('csi',csi,pos=1); assign('phi',phi,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); #################################### } ###### >>>>> End of CUBEFORSIM program <<<<< ######################## #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part VIII: Main CUBs functions #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ######################################## ### CUB00 ==> Main CUB(0,0) function ### ######################################## cub00<-function(ordinal,makeplot=TRUE){### Alternatively, makeplot=FALSE maxiter=500; toler=1e-6 ### Modify, if necessary tt0=proc.time(); serie=1:m; freq=tabulate(ordinal,nbins=m); n=sum(freq); aver=mean(ordinal); varcamp=mean(ordinal^2)-aver^2; ####################################################### ### ATTENTION !!!! this solution is not efficient !!! ### pai=0.5; csi=(m-aver)/(m-1); print(c(pai,csi)) ####################################################### inipaicsi=inibest(freq); pai=inipaicsi[1]; csi=inipaicsi[2]; ### pai=runif(1);csi=runif(1); Alternatively, random initial values ################################################################## loglik=loglikcub00(freq,pai,csi); # ******************************************************************** # ************* E-M algorithm for CUB(0,0) *************************** # ******************************************************************** nniter=1; while(nniter<=maxiter){ likold=loglik; bb=probbit(m,csi); aa=(1-pai)/(m*pai*bb); tau=1/(1+aa); ft=freq*tau; averpo=(t(serie)%*%ft)/sum(ft); pai=(t(freq)%*%tau)/n; # updated pai estimate csi=(m-averpo)/(m-1); # updated csi estimate if(csi<0.001){ csi=0.001;nniter=maxiter-1 } # print(c(pai,csi)); loglik=loglikcub00(freq,pai,csi); liknew=loglik; testll=abs(liknew-likold); ###### print(testll); # OPTIONAL printing: print(cbind(nniter,testll,pai,csi)); if(testll<=toler) break else {loglik=liknew}; # OPTIONAL printing: print(loglik); nniter=nniter+1; } ###### if(csi>0.999) csi=0.99; if(csi<0.001) csi=0.01; ### to avoid division by 0 !!! if(pai<0.001) pai=0.01; ### to avoid division by 0 !!! ### to ensure identifiability !!! ###### AICCUB00=-2*loglik+2*(2); BICCUB00=-2*loglik+log(n)*(2); ############################### varmat=varcovcub00(ordinal,pai,csi) ### Computation of var-covar of estimates ddd=diag(sqrt(1/diag(varmat))); #################################################################### # Print CUB(0,0) results of ML estimation #################################################################### nomi=rbind("pai","csi");stime=round(rbind(pai,csi),5); errstd=round(sqrt(diag(varmat)),5);wald=round(stime/errstd,5); pval=round(2*(1-pnorm(abs(wald))),20); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat("=====>>> C U B (0,0) model <<<===== ML-estimates via E-M algorithm ","\n"); cat("=======================================================================","\n"); cat("*** m=", m," *** Sample size: n=", n," *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:2){ cat(nomi[i]," ",stime[i]," ",errstd[i]," ",wald[i]," ",pval[i],"\n"); } cat("=======================================================================","\n"); cat("Variance-covariance matrix","\n"); print(round(varmat,5)); cat("=======================================================================","\n"); cat("Parameters correlation matrix","\n"); print(round(ddd%*%varmat%*%ddd,5)); ############################################################################## expcub=expcub00(m,pai,csi); varcub=varcub00(m,pai,csi); quota=(1-pai)/m; esq=sqrt(varmat[1,1])/m; lowq=quota-1.96*esq; uppq=quota+1.96*esq; llunif=-n*log(m); csisb=(m-aver)/(m-1); llsb=loglikcub00(tabulate(ordinal,nbins=m),1,csisb); epsilon=0.000001; logsat=-n*log(n)+sum((freq+epsilon)*log(freq+epsilon)); ### to avoid 0*log(0) devian=2*(logsat-loglik); cat("=======================================================================","\n"); cat("Log-lik(pai^,csi^) =",round(loglik,digits=8),"\n"); cat("Mean Log-likelihood=",round(loglik/n,digits=8),"\n"); cat("Log-lik(saturated) =",round(logsat,digits=8),"\n"); cat("Deviance =",round(devian,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("Log-lik(UNIFORM) =",round(llunif,digits=8),"\n"); cat("Log-lik(Shifted-BINOMIAL)=",round(llsb,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("AIC-CUB00 =",round(AICCUB00,digits=8),"\n"); cat("BIC-CUB00 =",round(BICCUB00,digits=8),"\n"); cat("=======================================================================","\n"); cat("Uncertainty Share: {(1-pai)/m} = ", quota,"(",esq,")","\n"); cat("Asympt.Conf.Int.(95%) =[",lowq,"; ",uppq,"]","\n"); cat("=======================================================================","\n"); theorpr=probcub00(m,pai,csi); pearson=((freq-n*theorpr))/sqrt(n*theorpr); X2=sum(pearson^2); relares=(freq/n-theorpr)/theorpr; LL2=1/(1+mean((freq/(n*theorpr)-1)^2)); II2=(loglik-llunif)/(logsat-llunif); diss00=0.5*sum(abs(theorpr-freq/n)); FF2=1-diss00; aupcub00=aup(m,freq,theorpr); mat1=cbind(pai,csi,loglik,n,X2,diss00); cat("Pearson Fitting measure ==> X^2 =",X2,"(p-val.=",1-pchisq(X2,m-3),")","\n"); cat("Lik-based fitting measure ==> L^2 =",LL2,"\n"); cat("Relative Log-lik index ==> I =",round(II2,digits=5),"\n"); cat("F^2 fitting measure ==> F^2 =",round(FF2,digits=5),"\n"); cat("Normed Dissimilarity index ==> Diss=",diss00,"\n"); cat("Area Under hit-Profile ==> AUP =",aupcub00,"\n"); cat("=======================================================================","\n"); cat("Observed average =", aver," Sample variance =",varcamp,"\n"); cat("Expectation of R~CUB(0,0) =", expcub," Variance of R~CUB(0,0) =",varcub,"\n"); cat("=======================================================================","\n"); cat("(R=r) Observed CUB-prob","Pearson","Relative res.","\n"); stampa=cbind(1:m,freq/n,theorpr,pearson,relares); print(stampa,digits=5); ################################################################ # Assignments as global variables: assign('name',value,pos=1);## ################################################################ assign('pai',pai,pos=1); assign('csi',csi,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); assign('diss',diss00,pos=1); assign('L2',LL2,pos=1); assign('II2',II2,pos=1); assign('n',n,pos=1); assign('AICCUB00',AICCUB00,pos=1); assign('BICCUB00',BICCUB00,pos=1); assign('aupcub00',aupcub00,pos=1); ### assign('nniter',nniter,pos=1); ################################################################ # /* CUB(0,0) ===>>> Only 1 plot with Diss=... */ ################################################################ rdiss00=round(diss00,3); ### Check on plotting ### if(makeplot==TRUE) DO PLOT; else (makeplot==FALSE) ===> NOPLOT ### ===>>> ## stringtitle=match.call()[[2]]; ### it writes "ordinal" if(makeplot==TRUE){ stringtitle="CUB model"; plot(cbind(1:m,1:m),cbind(theorpr,(freq/n)),las=1, main=paste(stringtitle, " (Diss =",round(diss00,digits=4),")"), xlim=c(1,m),ylim=c(0.0,1.1*max(theorpr,(freq/n))), xlab="Ordinal values of R=1,2,...,m", ylab=expression(paste("Observed relative frequencies (dots) and fitted probabilities (circles)"))); ### points(1:m,theorpr,pch=21,cex=1.5,lwd=2.0,type="b",lty=3); ### ex pch=8,col="red" points(1:m,freq/n,pch=16,cex=1.25,lwd=1.5); abline(h=0) } durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ###### >>>>> End of CUB(0,0) <<<<< ################################### ######################################## ### CUBp0 ==> Main CUB(p,0) function ### ######################################## cubp0<-function(ordinal,Y){ tt0=proc.time(); n=length(ordinal); p=ncol(as.matrix(Y)); aver=mean(ordinal); varcamp=mean(ordinal^2)-aver^2; YY=cbind(1,Y); # add 1's first column to Y ############################################################## ### previous >>> betjj=rep(0.1,p+1); csijj=(m-aver)/(m-1); have been modified ### >>> Now, it uses inibest(freq); ################################################################## serie=1:m; freq=tabulate(ordinal,nbins=m); inipaicsi=inibest(freq); betjj=rep(0.1,p+1); csijj=inipaicsi[2]; ############################################################## loglikjj=loglikcubp0(ordinal,Y,betjj,csijj); # ******************************************************************** # ************* E-M algorithm for CUB(p,0) *************************** # ******************************************************************** nniter=1; while(nniter<=maxiter){ loglikold=loglikjj; bb=probbit(m,csijj); vettn=bb[ordinal]; # probbit for all ordinal (r_i,i=1,2,...,n) aai=exp(-(YY%*%betjj)); ttau=1/(1+aai/(m*vettn)); # tau is a reserved word in R averpo=sum(ordinal*ttau)/sum(ttau); ################ effe1 is QPAI #################################### effe1<-function(bet,esterno=cbind(ttau,YY)){ tauno=esterno[,1]; covar=esterno[,2:ncol(esterno)]; sum(log(1+exp(-covar%*%bet))+(1-tauno)*(covar%*%bet)); } ################################## maximize w.r.t. bet ######## nlmaxbet=nlm(effe1,betjj); ################################################################ betjj=nlmaxbet$estimate; #updated bet estimates csijj=(m-averpo)/(m-1); #updated csi estimate loglikjj=loglikcubp0(ordinal,Y,betjj,csijj); # print(c(nniter,betjj,csijj,loglikjj)); #OPTIONAL PRINTING OF ITERATIONS testll=abs(loglikjj-loglikold); if(testll<=toler) break else {loglikold=loglikjj}; nniter=nniter+1;} bet=betjj; csi=csijj; loglik=loglikjj; #################################################################### AICCUBp0=-2*loglik+2*(p+2); BICCUBp0=-2*loglik+log(n)*(p+2); #################################################################### # Compute asymptotic standard errors of ML estimates #################################################################### varmat=varcovcubp0(ordinal,Y,bet,csi); if(det(varmat)<=0) stop("information matrix NOT positive definite"); ddd=diag(sqrt(1/diag(varmat))); cormat=(ddd%*%varmat)%*%ddd; #################################################################### # Print CUB(p,0) results of ML estimation #################################################################### nomi=c(paste("beta",0:(length(bet)-1),sep="_"),"csi "); stime=round(c(bet,csi),5); errstd=round(sqrt(diag(varmat)),5); wald=round(stime/errstd,5); pval=round(2*(1-pnorm(abs(wald))),20); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat("=====>>> C U B (p,0) model <<<===== ML-estimates via E-M algorithm ","\n"); cat("=======================================================================","\n"); cat(" Covariates for pai ==> p=", p,"\n"); cat("=======================================================================","\n"); cat("*** m=", m," *** Sample size: n=", n," *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:length(nomi)){ cat(nomi[i]," ",stime[i]," ",errstd[i]," ",wald[i]," ",pval[i],"\n"); } #################################################################### cat("=======================================================================","\n"); cat(" Variance-covariance matrix","\n"); rownames(varmat)=nomi;colnames(varmat)=nomi; print(round(varmat,5)); cat("=======================================================================","\n"); cat(" Parameters correlation matrix","\n"); rownames(cormat)=nomi;colnames(cormat)=nomi; print(round(cormat,5)); ############################################################################## cat("=======================================================================","\n"); cat("Log-lik(beta^,csi^) =",round(loglik,digits=8),"\n"); cat("Mean Log-likelihood =",round(loglik/n,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("AIC-CUBp0 =",round(AICCUBp0,digits=8),"\n"); cat("BIC-CUBp0 =",round(BICCUBp0,digits=8),"\n"); cat("=======================================================================","\n"); ################################################################ # Assignments as global variables ################################################################ assign('bet',bet,pos=1); assign('csi',csi,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); assign('n',n,pos=1); assign('AICCUBp0',AICCUBp0,pos=1); assign('BICCUBp0',BICCUBp0,pos=1); ### Comparing and plotting distributions if covariate for pai is dichotomus (0,1) ##### if(p==1 & length(unique(Y))==2) dicopai(ordinal,Y,m,bet,csi); ################################################################ durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ###### >>>>> End of CUB(p,0) <<<<< ################################### ######################################## ### CUB0q ==> Main CUB(0,q) function ### ######################################## cub0q<-function(ordinal,W){ tt0=proc.time(); n=length(ordinal); q=ncol(as.matrix(W)); aver=mean(ordinal); WW=cbind(1,W); # add 1's first column to W ############################################################## paijj=0.5; gamajj=inibestgama(ordinal,W); ### gamajj=rep(0.1,q+1); Alternatively !!! ############################################################## loglikjj=loglikcub0q(ordinal,W,paijj,gamajj); # ******************************************************************** # ************* E-M algorithm for CUB(0,q) *************************** # ******************************************************************** nniter=1; while(nniter<=maxiter){ loglikold=loglikjj; vettn=bitgama(ordinal,W,gamajj); ttau=1/(1+(1-paijj)/(m*paijj*vettn)); # tau is a reserved word in R ################ effe2 is QCSI #################################### effe2<-function(gama,esterno=cbind(ttau,ordinal,WW)){ tauno=esterno[,1]; ordd=esterno[,2]; covar=esterno[,3:ncol(esterno)]; sum(ttau*((ordinal-1)*(WW%*%gama)+(m-1)*log(1+exp(-WW%*%gama)))); } ################################# maximize w.r.t. gama ######## nlmaxgama=nlm(effe2,gamajj); ################################################################ gamajj=nlmaxgama$estimate; #updated gama estimates paijj=sum(ttau)/n; #updated pai estimate loglikjj=loglikcub0q(ordinal,W,paijj,gamajj); # print(c(nniter,paijj,gamajj,loglikjj)); #OPTIONAL PRINTING OF ITERATIONS testll=abs(loglikjj-loglikold); if(testll<=toler) break else {loglikold=loglikjj}; nniter=nniter+1;} pai=paijj; gama=gamajj; loglik=loglikjj; #################################################################### AICCUB0q=-2*loglik+2*(q+2); BICCUB0q=-2*loglik+log(n)*(q+2); #################################################################### # Compute asymptotic standard errors of ML estimates #################################################################### varmat=varcovcub0q(ordinal,W,pai,gama); if(det(varmat)<=0) stop("information matrix NOT positive definite"); ddd=diag(sqrt(1/diag(varmat))); cormat=(ddd%*%varmat)%*%ddd; #################################################################### ### Print CUB(0,q) results of ML estimation #################################################################### nomi=c("pai ",paste("gamma",0:(length(gama)-1),sep="_")); stime=round(c(pai,gama),5); errstd=round(sqrt(diag(varmat)),5); wald=round(stime/errstd,5); pval=round(2*(1-pnorm(abs(wald))),20); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat("=====>>> C U B (0,q) model <<<===== ML-estimates via E-M algorithm ","\n"); cat("=======================================================================","\n"); cat(" Covariates for csi ==> q=", q,"\n"); cat("=======================================================================","\n"); cat("*** m=", m," *** Sample size: n=", n," *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:length(nomi)){ cat(nomi[i]," ",stime[i]," ",errstd[i]," ",wald[i]," ",pval[i],"\n"); } #################################################################### cat("=======================================================================","\n"); cat(" Variance-covariance matrix","\n"); rownames(varmat)=nomi;colnames(varmat)=nomi; print(round(varmat,5)); cat("=======================================================================","\n"); cat(" Parameters correlation matrix","\n"); rownames(cormat)=nomi;colnames(cormat)=nomi; print(round(cormat,5)); ############################################################################## cat("=======================================================================","\n"); cat("Log-lik(pai^,gamma^) =",round(loglik,digits=8),"\n"); cat("Mean Log-likelihood =",round(loglik/n,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("AIC-CUB0q =",round(AICCUB0q,digits=8),"\n"); cat("BIC-CUB0q =",round(BICCUB0q,digits=8),"\n"); cat("=======================================================================","\n"); ################################################################ # Assignments as global variables ################################################################ assign('pai',pai,pos=1); assign('gama',gama,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); assign('n',n,pos=1); assign('AICCUB0q',AICCUB0q,pos=1); assign('BICCUB0q',BICCUB0q,pos=1); ### Comparing and plotting distributions if covariate for csi is dichotomus (0,1) ##### if(q==1 & length(unique(W))==2) dicocsi(ordinal,W,m,pai,gama); ################################################################ durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ###### >>>>> End of CUB(0,q) <<<<< ################################### ######################################## ### CUBpq ==> Main CUB(p,q) function ### ######################################## cubpq<-function(ordinal,Y,W){ tt0=proc.time(); n=length(ordinal); p=ncol(as.matrix(Y)); q=ncol(as.matrix(W)); aver=mean(ordinal); YY=cbind(1,Y); # add 1's first column to Y WW=cbind(1,W); # add 1's first column to W ################################################################################# ### previous >>> betjj=rep(0.1,p+1); gamajj=rep(0.1,q+1); have been modified ### >>> Now, it uses inibestgama(ordinal,W); ################################################################################# betjj=rep(0.1,p+1); gamajj=inibestgama(ordinal,W); ### Attention !!! ################################################################################# loglikjj=loglikcubpq(ordinal,Y,W,betjj,gamajj); # ******************************************************************** # ************* E-M algorithm for CUB(p,q) *************************** # ******************************************************************** nniter=1; while(nniter<=maxiter){ loglikold=loglikjj; vettn=bitgama(ordinal,W,gamajj); aai=exp(-(YY%*%betjj)); ttau=1/(1+aai/(m*vettn)); # tau is a reserved word in R ################ effe1 is QPAI #################################### effe1<-function(bet,esterno=cbind(ttau,YY)){ tauno=esterno[,1]; covar=esterno[,2:ncol(esterno)]; sum(log(1+exp(-covar%*%bet))+(1-tauno)*(covar%*%bet)); } ################ effe2 is QCSI #################################### effe2<-function(gama,esterno=cbind(ttau,ordinal,WW)){ tauno=esterno[,1]; ordd=esterno[,2]; covar=esterno[,3:ncol(esterno)]; sum(ttau*((ordinal-1)*(WW%*%gama)+(m-1)*log(1+exp(-WW%*%gama)))); } #################### maximize w.r.t. bet and gama ############ nlmaxbet=nlm(effe1,betjj); nlmaxgama=nlm(effe2,gamajj); ################################################################ betjj=nlmaxbet$estimate; #updated bet estimates gamajj=nlmaxgama$estimate; #updated gama estimates loglikjj=loglikcubpq(ordinal,Y,W,betjj,gamajj); # print(c(nniter,betjj,gamajj,loglikjj)); #OPTIONAL PRINTING OF ITERATIONS testll=abs(loglikjj-loglikold); if(testll<=toler) break else {loglikold=loglikjj}; nniter=nniter+1;} bet=betjj; gama=gamajj; loglik=loglikjj; #################################################################### AICCUBpq=-2*loglik+2*(p+q+2); BICCUBpq=-2*loglik+log(n)*(p+q+2); #################################################################### # Compute asymptotic standard errors of ML estimates #################################################################### varmat=varcovcubpq(ordinal,Y,W,bet,gama); if(det(varmat)<=0) stop("information matrix NOT positive definite"); ddd=diag(sqrt(1/diag(varmat))); cormat=(ddd%*%varmat)%*%ddd; #################################################################### # Print CUB(p,q) results of ML estimation #################################################################### nomi=c(paste("beta",0:(length(bet)-1),sep="_"),paste("gamma",0:(length(gama)-1),sep="_")); stime=round(c(bet,gama),5); errstd=round(sqrt(diag(varmat)),5); wald=round(stime/errstd,5); pval=round(2*(1-pnorm(abs(wald))),20); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat("=====>>> C U B (p,q) model <<<===== ML-estimates via E-M algorithm ","\n"); cat("=======================================================================","\n"); cat(" Covariates for pai ==> p=",p," and Covariates for csi ==> q=", q,"\n"); cat("=======================================================================","\n"); cat("*** m=", m," *** Sample size: n=", n," *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:length(nomi)){ cat(nomi[i]," ",stime[i]," ",errstd[i]," ",wald[i]," ",pval[i],"\n"); } #################################################################### cat("=======================================================================","\n"); cat(" Variance-covariance matrix","\n"); rownames(varmat)=nomi;colnames(varmat)=nomi; print(round(varmat,5)); cat("=======================================================================","\n"); cat(" Parameters correlation matrix","\n"); rownames(cormat)=nomi;colnames(cormat)=nomi; print(round(cormat,5)); ############################################################################## cat("=======================================================================","\n"); cat("Log-lik(beta^,gamma^)=",round(loglik,digits=8),"\n"); cat("Mean Log-likelihood =",round(loglik/n,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("AIC-CUBpq =",round(AICCUBpq,digits=8),"\n"); cat("BIC-CUBpq =",round(BICCUBpq,digits=8),"\n"); cat("=======================================================================","\n"); ### Comparing and plotting distributions if the same covariate for pai and csi ### is dichotomus (0,1) ##### if(p==1 & length(unique(Y))==2 & q==1 & length(unique(W))==2) dicopaicsi(ordinal,Y,W,m,bet,gama); ################################################################ # Assignments as global variables ################################################################ assign('bet',bet,pos=1); assign('gama',gama,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); assign('n',n,pos=1); assign('AICCUBpq',AICCUBpq,pos=1); assign('BICCUBpq',BICCUBpq,pos=1); ################################################################ durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ###### >>>>> End of CUB(p,q) <<<<< ################################ ####################################### ### CUBSHE ==> Main CUBSHE function ### ####################################### cubshe<-function(ordinal,ccc){ tt0=proc.time(); serie=1:m; freq=tabulate(ordinal,nbins=m); n=sum(freq); aver=mean(ordinal); varcamp=mean(ordinal^2)-aver^2; dd=rep(0,m);dd[ccc]=1; ### Alternatively: dd=ifelse(serie==ccc,1,0) ##################################################################### # ATTENTION !!!! Estimation routines use the function 'probcubshe1' ##################################################################### ################################################################ # !!! Initial values; we are improving them !!! ################################################################ ### pai1=1/3; pai2=1/3; csi=(m-aver)/(m-1);### implies delta=1/3 pai1=0.45; pai2=0.45; csi=(m-aver)/(m-1); ### implies delta=0.1 ################################################################ loglik=loglikcubshe(freq,pai1,pai2,csi,ccc); # ******************************************************************** # ************* E-M algorithm for CUBSHE ***************************** # ******************************************************************** nniter=1; while(nniter<=maxiter){ likold=loglik; bb=probbit(m,csi); tau1=pai1*bb; tau2=pai2*(1/m); denom=tau1+tau2+(1-pai1-pai2)*dd; tau1=tau1/denom; tau2=tau2/denom; tau3=1-tau1-tau2; numaver=sum(serie*freq*tau1); denaver=sum(freq*tau1); averpo=numaver/denaver; pai1=sum(freq*tau1)/n; #updated pai1 estimate pai2=sum(freq*tau2)/n; #updated pai2 estimate csi=(m-averpo)/(m-1); #updated csi estimate if(csi<0.001){ csi=0.001;nniter=maxiter-1 } loglik=loglikcubshe(freq,pai1,pai2,csi,ccc); liknew=loglik; testll=abs(liknew-likold); #print(cbind(nniter,testll,pai1,pai2,csi,loglik)); if(testll<=toler) break else {loglik=liknew}; nniter=nniter+1; } ###### if(csi>0.999) csi=0.98; ### to avoid division by 0 !!! if(csi<0.001) csi=0.02; ### to avoid division by 0 !!! if(pai1<0.001) pai1=0.02; ### to ensure identifiability !!! #################################################################### # Compute asymptotic standard errors of ML estimates #################################################################### varmat=varcovcubshe(pai1,pai2,csi,ccc,n); ddd=diag(sqrt(1/diag(varmat))); ###### delta=1-pai1-pai2; esdelta=sqrt(varmat[1,1]+varmat[2,2]+2*varmat[1,2]); pvaldelta=round(2*(1-pnorm(abs(delta/esdelta))),20); paistar=pai1/(pai1+pai2); espaistar=sqrt((pai1^2*varmat[2,2]+pai2^2*varmat[1,1]-2*pai1*pai2*varmat[1,2]))/(pai1+pai2)^2; pvalpaistar=round(2*(1-pnorm(abs(paistar/espaistar))),20); #################################################################### # Print CUBSHE results of ML estimation #################################################################### nomi=rbind("pai1","pai2","csi"); stime=rbind(pai1,pai2,csi); errstd=sqrt(diag(varmat)); wald=stime/errstd; pval=round(2*(1-pnorm(abs(wald))),20); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat("==>>> C U B + SHELTER model <<<===== ML-estimates via E-M algorithm ","\n"); cat("=======================================================================","\n"); cat(" m=", m," shelter=", ccc," Sample size: n=", n, " *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:3){ cat(nomi[i]," ",round(stime[i],5)," ",round(errstd[i],5)," ",round(wald[i],5)," ",pval[i],"\n"); } cat("=======================================================================","\n"); cat("delta"," ",delta, " ",esdelta, " ",delta/esdelta," ",pvaldelta,"\n"); cat("(1-pai1-pai2)","\n"); cat("=======================================================================","\n"); cat("paistar"," ",paistar," ",espaistar," ",paistar/espaistar," ",pvalpaistar,"\n"); cat("[pai1/(pai1+pai2)]","\n"); cat("=======================================================================","\n"); cat("Variance-covariance matrix","\n"); print(round(varmat,5)); cat("=======================================================================","\n"); cat("Parameters correlation matrix","\n"); print(round(ddd%*%varmat%*%ddd,5)); llunif=-n*log(m); csisb=(m-aver)/(m-1); llsb=loglikcub00(tabulate(ordinal,nbins=m),1,csisb); cat("=======================================================================","\n"); cat("Log-lik(theta^) =",round(loglik,digits=8),"\n"); cat("Mean Log-likelihood=",round(loglik/n,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("Log-lik(UNIFORM) =",round(llunif,digits=8),"\n"); cat("Log-lik(Shifted-BINOMIAL)=",round(llsb,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); theorpr=probcubshe1(m,pai1,pai2,csi,ccc); dissshe=dissom(theorpr,freq/n); llunif=-n*log(m); epsilon=0.000001; logsat=-n*log(n)+sum((freq+epsilon)*log(freq+epsilon)); ### to avoid 0*log(0) X2=sum(((n*theorpr-freq)^2)/(n*theorpr)); LL2=1/(1+mean((freq/(n*theorpr)-1)^2)); II2=(loglik-llunif)/(logsat-llunif); FF2=1-dissshe; aupcubshe=aup(m,freq,theorpr); cat("Pearson Fitting measure ==> X^2 =",X2,"(p-val.=",1-pchisq(X2,m-3),")","\n"); cat("Normed Dissimilarity index ==> Diss=",round(dissshe,digits=5),"\n"); cat("F^2 fitting measure ==> F^2 =",round(FF2,digits=5),"\n"); cat("Lik-based fitting measure ==> L^2 =",LL2,"\n"); cat("Relative Log-lik index ==> I =",round(II2,digits=5),"\n"); cat("Area Under hit-Profile ==> AUP =",round(aupcubshe,digits=5),"\n"); cat("=======================================================================","\n"); ################################################################ # Assignments as global variables ################################################################ assign('pai1',pai1,pos=1); assign('pai2',pai2,pos=1); assign('csi',csi,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); assign('dissshw',dissshe,pos=1); assign('n',n,pos=1); assign('delta',delta,pos=1); assign('esdelta',esdelta,pos=1); assign('paistar',paistar,pos=1); ################################################################ # /* CUB-she ===>>> Only 1 plot with Diss=... */ ################################################################ plot(cbind(1:m,1:m),cbind(theorpr,(freq/n)), main=paste("CUB model with shelter effect"," (Diss =",round(dissshe,digits=4),")"), xlim=c(1,m),ylim=c(0,1.1*max(theorpr,(freq/n))),las=1, xlab="Ordinal values of r=1,2,...,m", ylab=expression(paste("Observed relative frequencies (dots) and fitted probabilities (circles)"))); points(1:m,theorpr,pch=21,cex=1.5); points(1:m,freq/n,pch=16,cex=1.2); ### points(ccc,theorpr[ccc]-delta,pch=8); abline(h=0); durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ###### >>>>> End of CUBSHE(0,0) <<<<< ################################### ############################################################## ############################################################## ### >>>>>>>>>>>> Main CUBE000 function <<<<<<<<<<<<<<<<<<< ### ############################################################## cube000<-function(ordinal,starting=rep(0.1,3),maxiter=500, toler=1e-6,makeplot=TRUE,expinform=FALSE){ tt0=proc.time(); freq=tabulate(ordinal,nbins=m); n=sum(freq); aver=mean(ordinal); varcamp=mean(ordinal^2)-aver^2; ######################################################## #* Defaults: toler = 1e-06; maxiter = 500; ***** ######################################################## #(0)# initial estimates, not efficient: ### pai=starting[1]; csi=starting[2]; phi=starting[3]; ### Alternatively, more efficient: starting=inibestcube(ordinal) pai=starting[1]; csi=starting[2]; phi=starting[3]; #(0)# log-lik loglik=loglikcube(m,pai,csi,phi,freq); # ******************************************************************** # ************* E-M algorithm for CUBE ************************* # ******************************************************************** nniter=1; while(nniter<=maxiter){ likold=loglik; #(1)# betar bb=betar(m,csi,phi); aa=(1-pai)/(m*pai*bb); #(2)# taunor tauno=1/(1+aa); #(3)# pai(k+1) pai=sum(freq*tauno)/n; # updated pai estimate paravecjj=c(csi,phi); #(4)# Q(k+1) dati=cbind(tauno,freq); ################ EFFECUBE is Q(csi,phi) ########################### #(5)# (csi(k+1),phi(k+1)) ################################## maximize w.r.t. paravec ######## optimestim=optim(paravecjj,EFFECUBE,dati=dati); # print(nlmaxg) ################################################################ #(6)# theta(k+1) paravecjj=optimestim$par; # updated paravec estimates csi=paravecjj[1]; phi=paravecjj[2]; ########################################## # if(pai<0.001){pai=0.001; nniter=maxiter-1} # if(csi<0.001){csi=0.001; nniter=maxiter-1} # if(phi<0.001){phi=0.001; nniter=maxiter-1} # if(pai>0.999){pai=0.98} ### to avoid division by 0 !!! # if(csi>0.999){csi=0.98} ### to avoid division by 0 !!! ###################################### print(c(nniter,pai,csi,phi)); #(7)# elle(theta(k+1)) liknew=loglikcube(m,pai,csi,phi,freq); #(8)# test testll=abs(liknew-likold); # OPTIONAL printing: print(testll); # OPTIONAL printing: print(cbind(nniter,testll,pai,csi,phi)); if(testll<=toler) break else {loglik=liknew}; # OPTIONAL printing: print(loglik); nniter=nniter+1; } loglik=liknew; ###### End of E-M algorithm for CUBE *********************************************** ###### AICCUBE=-2*loglik+2*(3); BICCUBE=-2*loglik+log(n)*(3); # ******************************************************** # Compute ML var-cov matrix and print result for CUBE # ******************************************************** if(expinform==TRUE){ varmat=varcovcubeexp(m,pai,csi,phi,n) } else{ varmat=varcovcubeobs(m,pai,csi,phi,freq) } #################################################################### # Print CUBE results of ML estimation #################################################################### nomi=rbind("pai ","csi ","phi "); stime=rbind(pai,csi,phi); errstd=sqrt(diag(varmat)); wald=stime/errstd; pval=round(2*(1-pnorm(abs(wald))),20); ddd=diag(sqrt(1/diag(varmat))); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat("> CUBE model Inference ","\n"); cat("=======================================================================","\n"); if(expinform==TRUE){ cat("Maximum Likelihood estimates (E-M algorithm) ..... Expected Information","\n"); } else{ cat("Maximum Likelihood estimates (E-M algorithm) ..... Observed Information","\n"); } cat("=======================================================================","\n"); cat("*** m=", m," *** Sample size: n=", n," *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:3){ cat(nomi[i]," ",round(stime[i],5)," ",round(errstd[i],5)," ",round(wald[i],5)," ",pval[i],"\n"); } cat("=======================================================================","\n"); cat("Variance-covariance matrix","\n"); print(round(varmat,5)); cat("=======================================================================","\n"); cat("Parameters correlation matrix","\n"); print(round(ddd%*%varmat%*%ddd,5)); ### Log-likelihood comparisons ############################################## llunif=-n*log(m); csisb=(m-aver)/(m-1); llsb=loglikcub00(freq,1,csisb); epsilon=0.000001; logsat=-n*log(n)+sum((freq+epsilon)*log(freq+epsilon)); ### to avoid 0*log(0) devian=2*(logsat-loglik); cat("\n"); cat("==================================================== Log-likelihoods ===","\n"); cat("Log-lik(theta^) =",round(loglik,digits=8),"\n"); cat("Log-lik(Saturated) =",round(logsat,digits=8),"\n"); cat("Log-lik(Shifted-Binomial)=",round(llsb,digits=8),"\n"); cat("Log-lik(Uniform) =",round(llunif,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("Mean Log-likelihood=",round(loglik/n,digits=8),"\n"); cat("Deviance =",round(devian,digits=8),"\n"); ### Fitting measures ######################################################### cat("=========================================== Global fitting measures ===","\n"); cat("\n"); theorpr=probcube(m,pai,csi,phi); dissom=0.5*sum(abs(theorpr-freq/n)); aupcube=aup(m,freq,theorpr); pearson=((freq-n*theorpr))/sqrt(n*theorpr); X2=sum(pearson^2); relares=(freq/n-theorpr)/theorpr FF2=1-dissom; LL2=1/(1+mean((freq/(n*theorpr)-1)^2)); II2=(loglik-llunif)/(logsat-llunif); cat("Pearson Fitting measure ==> X^2 =",X2,"(p-val.=",1-pchisq(X2,m-4),")","\n"); cat("F^2 fitting measure ==> F^2 =",round(FF2,digits=5),"\n"); cat("Normed Dissimilarity index ==> Diss=",round(dissom,digits=5),"\n"); cat("Area Under hit-Profile ==> AUP =",aupcube,"\n"); cat("Lik-based fitting measure ==> L^2 =",round(LL2,digits=5),"\n"); cat("Relative Log-lik index ==> I =",round(II2,digits=5),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("AIC-CUBE =",round(AICCUBE,digits=8),"\n"); cat("BIC-CUBE =",round(BICCUBE,digits=8),"\n"); cat("=======================================================================","\n"); cat("(R=r) Observed CUBE-prob","Pearson","Relative res.","\n"); stampa=cbind(1:m,freq/n,theorpr,pearson,relares); print(stampa,digits=5); ############################################################################# # /* CUBE ===> plot with Diss=.......... optioned by makeplot=TRUE/FALSE */ ############################################################################# if(makeplot==TRUE){ stringtitle="CUBE model estimation "; plot(cbind(1:m,1:m),cbind(theorpr,(freq/n)),las=1, main=paste(stringtitle, " (Diss =",round(dissom,digits=4),")"), xlim=c(1,m),ylim=c(0.0,1.1*max(theorpr,(freq/n))), xlab="Ordinal values of R=1,2,...,m", ylab="Observed relative frequencies (dots) and fitted probabilities (circles)"); ### points(1:m,theorpr,pch=21,cex=1.5,lwd=2.0,type="b",lty=3); ### ex pch=8,col="red" points(1:m,freq/n,pch=16,cex=1.25,lwd=1.5); abline(h=0) } ##################################################################### # Assignments as global variables: assign('name',value,pos=1); ##################################################################### assign('nniter',nniter,pos=1); assign('pai',pai,pos=1); assign('csi',csi,pos=1); assign('phi',phi,pos=1); assign('varmat',varmat,pos=1); assign('loglik',loglik,pos=1); assign('BICCUBE',BICCUBE,pos=1); assign('FF2',FF2,pos=1); assign('II2',II2,pos=1); assign('aupcube',aupcube,pos=1); assign('n',n,pos=1); #################################### durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ###### >>>>> End of CUBE program <<<<< ################################### ############################################################## #>>>>>>>>>>>>>>> Main CUBECOV function <<<<<<<<<<<<<<<<<<<<<<< ############################################################## cubecov<-function(ordinal,Y,W,Z,starting=rep(0.1,ncol(cbind(Y,W,Z)+3)), maxiter=500,toler=1e-6){ tt0=proc.time(); n=length(ordinal); p=ncol(as.matrix(Y)); q=ncol(as.matrix(W)); v=ncol(as.matrix(Z)); aver=mean(ordinal); YY=cbind(1,Y); # add 1's first column to Y WW=cbind(1,W); # add 1's first column to W ZZ=cbind(1,Z); # add 1's first column to Z # ********************************** # *** E-M algorithm for CUBECOV *** # ********************************** ################################################################# # 00.# Initial values of parameters ............. Attention !!! # ################################################################# betjj=starting[1:(p+1)]; gamajj=starting[(p+2):(p+q+2)]; alphajj=starting[(p+q+3):(p+q+v+3)]; ################################ # * * * Iterative Loop * * *# ################################ nniter=1; while(nniter<=maxiter){ ################################################################# # 01.# Initial values of vectors i=1,2,..,n ################################################################# paijj=logis(Y,betjj); csijj=logis(W,gamajj); phijj=logis(Z,alphajj); ################################################################# # 02.# Computation of beta-binomial distribution i=1,2,..,n ################################################################# betabin=betabinomial(m,csijj,phijj,ordinal); ################################################################# # 03.# Computation of CUBE probability distribution i=1,2,..,n ################################################################# probi=paijj*(betabin-1/m)+1/m likold=sum(log(probi)) ### loglikcubecov(ordinal,Y,W,Z,bet,gama,alpha)### log-likelihood ################################################################# # 4.# Computation of conditional probability i=1,2,..,n ################################################################# taui=1-(1-paijj)/(m*probi) ################################################################# # 5a.# Computation of Q_1(beta) ################################################################# ################ Q_1(beta) ######################### Quno<-function(bet,taui=taui,YY=YY){ tauno=taui; covar=YY; ybeta=covar%*%bet; -sum(tauno*ybeta-log(1+exp(ybeta))); ### change sign for optim } ################################################################# # 5b.# Computation of Q_2(gamma,alpha) ################################################################# ################# unify parameter vectors param=c(gamajj,alphajj) ################ Q_2(gamma,alpha) ################## Qdue<-function(param,taui=taui,ordinal=ordinal,W=W,Z=Z){ tauno=taui; ordin=ordinal; gama=param[1:(q+1)]; alpha=param[(q+2):(q+v+2)]; csi=logis(W,gama); phi=logis(Z,alpha); betabin=betabinomial(m,csi,phi,ordinal); -sum(tauno*log(betabin)); ### change sign for optim } ################################################################# # 6.# Maximization of Q_1(beta) and Q_2(beta) ################################################################# ### maximize w.r.t. bet and gama ######### optimbet =optim(betjj,Quno,taui=taui,YY=YY); optimparam=optim(param,Qdue,taui=taui,ordinal=ordinal,W=W,Z=Z); ################################################################# # 7.# Computation of updated estimates and log-likelihood ################################################################# betjj=optimbet$par; #updated bet estimates paramjj=optimparam$par; gamajj=paramjj[1:(q+1)]; #updated gama estimates alphajj=paramjj[(q+2):(q+v+2)]; #updated alpha estimates ### updated log-likelihood liknew=loglikcubecov(ordinal,Y,W,Z,betjj,gamajj,alphajj) ################################################################# # 8.# Checking improvement of updated log-likelihood ################################################################# # print(c(nniter,betjj,gamajj,loglikjj)); #OPTIONAL PRINTING OF ITERATIONS testloglik=abs(liknew-likold); print(nniter); print(round(c(liknew,likold,testloglik),digits=7)) if(testloglik<=toler) break else {likold=liknew}; nniter=nniter+1 } ################################################################# # 8.# Final ML estimates and maximized log-likelihood ################################################################# bet=betjj; gama=gamajj; alpha=alphajj; loglik=liknew; ###### End of E-M algorithm for CUBE *********************************************** paramest=c(bet,gama,alpha) # intermediate PRINTING # stampa=c(nniter,paramest,loglik); # print(round(stampa,digits=5)) ###### AICCUBE=-2*loglik+2*(p+q+v+3); BICCUBE=-2*loglik+log(n)*(p+q+v+3); ############################################################ # Compute asymptotic standard errors of ML CUBE estimates ## ############################################################ varmat=varcovcubecov(m,ordinal,Y,W,Z,bet,gama,alpha) if(det(varmat)<=0) stop("information matrix NOT positive definite"); ddd=diag(sqrt(1/diag(varmat))); cormat=(ddd%*%varmat)%*%ddd; ################################################## # Print CUBE-covariates results of ML estimation # ################################################## nomi=c(paste("beta",0:(length(bet)-1),sep="_"), paste("gamma",0:(length(gama)-1),sep="_"), paste("alpha",0:(length(alpha)-1),sep="_")); stime=round(c(bet,gama,alpha),5); errstd=round(sqrt(diag(varmat)),5); wald=round(stime/errstd,5); pval=round(2*(1-pnorm(abs(wald))),20); cat("\n"); cat("=======================================================================","\n"); cat("============== CUB Program: version 3.0 (September 2013) ==============","\n"); cat("=======================================================================","\n"); cat(">> C U B E model with covariates << ML-estimates via E-M algorithm ","\n"); cat("=======================================================================","\n"); cat("Covar.pai: p=", p,"; Covar.csi: q=", q,"; Covar.phi: v=",v, "\n"); cat("=======================================================================","\n"); cat("*** m=", m," *** Sample size: n=", n," *** Iterations=",nniter," (maxiter=500)","\n"); cat("=======================================================================","\n"); cat("parameters ML-estimates stand.errors Wald-test p-value ","\n"); cat("=======================================================================","\n"); for(i in 1:length(nomi)){ cat(nomi[i]," ",stime[i]," ",errstd[i]," ",wald[i]," ",pval[i],"\n"); } #################################################################### cat("=======================================================================","\n"); cat(" Variance-covariance matrix","\n"); rownames(varmat)=nomi;colnames(varmat)=nomi; print(round(varmat,5)); cat("=======================================================================","\n"); cat(" Parameters correlation matrix","\n"); rownames(cormat)=nomi;colnames(cormat)=nomi; print(round(cormat,5)); ############################################################################## cat("=======================================================================","\n"); cat("Log-lik(beta^,gamma^)=",round(loglik,digits=8),"\n"); cat("Mean Log-likelihood =",round(loglik/n,digits=8),"\n"); cat("-----------------------------------------------------------------------","\n"); cat("AIC-CUBE =",round(AICCUBE,digits=8),"\n"); cat("BIC-CUBE =",round(BICCUBE,digits=8),"\n"); cat("=======================================================================","\n"); ##################################################################### # Assignments as global variables: assign('name',value,pos=1); ##################################################################### assign('nniter',nniter,pos=1); assign('bet',bet,pos=1); assign('gama',gama,pos=1); assign('alpha',alpha,pos=1); assign('loglik',loglik,pos=1); assign('BICCUBE',BICCUBE,pos=1); assign('n',n,pos=1); assign('varmat',varmat,pos=1); #################################### durata=proc.time()-tt0;durata=durata[1]; cat("=======================================================================","\n"); cat("Elapsed time=",durata,"seconds","===============>>>",date(),"\n"); cat("=======================================================================","\n","\n"); } ### >>>>> End of CUBECOV function <<<<< ### #============================================================================== #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #** Part IX: Main call for CUB and CUBE models ........................ #** 같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같같 #============================================================================== ### ************************** Main call to CUB() ***************************** ################################################################################ #**............................................................................# #**............................................................................# #**............................................................................# #**........... ........................# #**........... ������� �� �� �����\ ........................# #**........... �� �� �� �� � .......=============....# #**........... �� �� �� �� � ....... Version 3.0 ....# #**........... �� �� �� ������� .......=============....# #**........... �� �� �� �� � ........................# #**........... �� �� �� �� � ........................# #**........... ������� �������� �����/ ........................# #**........... ........................# #**............................................................................# #**............................................................................# #**............................................................................# #** This function calla the main procedures for building statistical models # #** of the family of CUB models and extensions............................ # #**............................................................................# ################################################################################ CUB<-function(ordinal,Y=0,W=0,shelter=0){ maxiter=500; toler=0.000001; assign('maxiter',maxiter,pos=1); assign('toler',toler,pos=1); ry=nrow(as.matrix(Y)); rw=nrow(as.matrix(W)); shelter=as.numeric(shelter); if(shelter!=0) cubshe(ordinal,shelter) # Assigned variables are: pai1,pai2,csi,notpd,varmat,loglik,diss,n,delta,esdelta,paistar else{ if(ry==1 & rw==1) cub00(ordinal) # Assigned variables are: pai,csi,varmat,loglik,n,AICCUB00,BICCUB00 else{ if(ry!=1 & rw==1) cubp0(ordinal,Y) # Assigned variables are: bet,csi,varmat,loglik,n,AICCUBp0,BICCUBp0 else{ if(ry==1 & rw!=1) cub0q(ordinal,W) # Assigned variables are: pai,gama,varmat,loglik,n,AICCUB0q,BICCUB0q else{ if(ry!=1 & rw!=1) cubpq(ordinal,Y,W) # Assigned variables are: bet,gama,varmat,loglik,n,AICCUBpq,BICCUBpq else stop("wrong variables specification") } } } } } ################################################################################ ### ************************** Main call to CUBE() ***************************** ################################################################################ #**............................................................................# #**............................................................................# #**............................................................................# #**...... .......................# #**...... ������� �� �� �����\ ������� .......................# #**...... �� �� �� �� � � ......=============....# #**...... �� �� �� �� � � ...... Version 3.0 ....# #**...... �� �� �� ������� ������� ......=============....# #**...... �� �� �� �� � � .......................# #**...... �� �� �� �� � � .......................# #**...... ������� �������� �����/ ������� .......................# #**....... .......................# #**............................................................................# #**............................................................................# #**............................................................................# #** This function calls the main procedures for building CUBE models # #** without and with covariates ..............................................# #**............................................................................# ################################################################################ CUBE<-function(ordinal,Y=0,W=0,Z=0,starting=rep(0.1,20),maxiter=500,toler=1e-6, makeplot=FALSE,expinform=FALSE){ ry=nrow(as.matrix(Y)); rw=nrow(as.matrix(W)); rz=nrow(as.matrix(Z)); if(ry==1 & rw==1 & rz==1) {cube000(ordinal,starting=rep(0.1,3),maxiter,toler,makeplot=TRUE,expinform)} else{ cubecov(ordinal,Y,W,Z,starting=rep(0.1,ncol(cbind(Y,W,Z))+3),maxiter,toler) } } ################################################################################