/* Program for fitting a bivariate Plackett distribution with CUB margins to the olive oil dataset as in the article: M.CORDUAS (2014): ANALYZING BIVARIATE ORDINAL DATA WITH CUB MARGINS. INPUT FILE: OLIVE_OIL.DAT This is a GAUSS dataset containing a data matrix with 1000 rows and 3 columns. The data originates from 1000 interviews. It includes the variables: Colour, Flavour, Duminfo. "Colour" refers to rates expressed by interviewees about the colour of extra-virgin olive oil using a 7 points Likert scale. "Flavour" refers to the rates about extra-virgin oil flavour. "Duminfo" refers to the degree of knowledge of the consumer about the olive oil (low-medium or high level). ________________________________________________________________________________________________________________________________________ This program estimates the Plackett copula with CUB margins and applies it to the variables: (Colour, Flavour) from olive oil data set Colour follows a CUB model with the covariate Duminfo acting on the feeling parameter Flavour follows a CUB model with the covariate Duminfo on both the feeling and uncertainty parameter The association parameter of the Plackett copula is affected by Duminfo. ________________________________________________________________________________________________________________________________________ LIST of PROCS - freqdoppia computes a joint bivariate distribution - solplack computes the Plackett copula cumulative value - Likpsicov computes the log-likelihood of the model including one covariate - Hpro computes the joint cumulative distribution of the Plackett copula given the margins and the association parameter, and the Plackett probability distribution -OPENVARS open a gauss dataset The program uses the CML (constrained maximum likelihood module) developed by Ronald Schoenberg for Aptech GAUSS (see reference Manual). ______________________________________________________________________________ OUTPUT: the program generates the file: OUTRESGATE.OUT ____________________________________________________________________________ */ proc(1)=freqdoppia(m1,m2); local freqdop,a,i,m; m=varget("m"); freqdop=zeros(m,m); i=1; do while i<=m; a=selif(m1,m2.==i); freqdop[i,.]=counts(a,seqa(1,1,m))'/rows(m1); i=i+1; endo; retp(freqdop); endp; proc(1)=solplack(F,G,PSI); local S,Hpla,a; a=PSI-1; S=1+a*(F+G); Hpla=(S-sqrt(S^2 -4*psi*a*F*G))/(2*a); retp(Hpla); endp; clear m,F0,G0,F1,G1,likePAE; proc(1)=LIKPSICOV(para,fdopempi); local m,lik,lik0,como,HHMAT0,PP0,FF0,GG0,lik1,HHMAT1,PP1,FF1,GG1,psi0,psi1; m=varget("m"); psi0=exp(para[1]); psi1=exp(para[1]+para[2]); FF0=varget("F0"); GG0=varget("G0"); {HHMAT0,PP0}=Hpro(FF0,GG0,psi0); lik0=sumc(sumc(fdopempi[.,1:m].*ln(PP0))); FF1=varget("F1"); GG1=varget("G1"); {HHMAT1,PP1}=Hpro(FF1,GG1,psi1); lik1=sumc(sumc(fdopempi[.,m+1:2*m].*ln(PP1))); lik=lik1+lik0; como=varput(lik,"likePAE"); retp(lik); endp; proc(2)=Hpro(FF,GG,psi); local Hmat, j,h, p,i,m; m=rows(FF); Hmat=zeros(m,m); j=1; do while j<=m; h=1; do while h<=m; Hmat[j,h]=solplack(FF[j],GG[h],psi); h=h+1; endo; j=j+1; endo; p=zeros(m,m); p[1,1]=hmat[1,1]; p[1,2:m]=hmat[1,2:m]-hmat[1,1:m-1]; p[2:m,1]=hmat[2:m,1]-hmat[1:m-1,1]; i=2; do while i<=m; h=2; do while h<=m; p[i,h]=hmat[i,h]-sumc(sumc(p[1:i,1:h])); h=h+1; endo; i=i+1; endo; retp(Hmat,p); endp; proc(0)= OPENVARS(dataset); local xxx,varnames,yyy; xxx=loadd(dataset); varnames=getname(dataset); yyy=VARPUT(varnames,"varnames"); MAKEVARS(xxx,0,varnames); endp; clear Colour,Flavour,Duminfo,m1,m2,pa,likepae; call openvars("OLIVE_OIL"); dati=Colour~Flavour~Duminfo; output file=outresgate.out reset; dati0=SELIF(dati,DATI[.,3].==0); dati1=SELIF(dati,DATI[.,3].==1); __numplot=3; __m=7; m=7; mteo0=zeros(m,2); mteo1=zeros(m,2); mfre0=zeros(m,2); mfre1=zeros(m,2); med0=zeros(2,1); med1=zeros(2,1); medemp0=zeros(2,1); medemp1=zeros(2,1); csi0=zeros(2,1); pigre0=zeros(2,1); csi1=zeros(2,1); pigre1=zeros(2,1); pigre=zeros(2,1); psi=zeros(2,1); fatt1=DATI[.,3]; @Variable1@ mm=dati[.,1]; dd=saved(mm,"comodo1","m1"); call openvars("comodo1"); @::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: The program recall the CUB estimation procedure implemented in GAUSS and R according to the algorithm published in the article: Piccolo D."Observed information matrix for MUB models", QDS, 2006. This is available at https://www.researchgate.net/profile/Domenico_Piccolo2. :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::@ {mat1, mat2} = CUB("comodo1","m1",0,fatt1); csi0[1]=1/(1+exp(-mat1[2])); csi1[1]=1/(1+exp(-mat1[2]-mat1[3])); pigre[1]=mat1[1]; mteo0[.,1]=PROBCUB00(m,pigre[1],csi0[1]); mteo1[.,1]=PROBCUB00(m,pigre[1],csi1[1]); mfre0[.,1]=counts(dati0[.,1],seqa(1,1,m))/rows(dati0); mfre1[.,1]=counts(dati1[.,1],seqa(1,1,m))/rows(dati1); med0[1]=sumc(seqa(1,1,m).*mteo0[.,1]); med1[1]=sumc(seqa(1,1,m).*mteo1[.,1]); medemp0[1]=sumc(seqa(1,1,m).*mfre0[.,1]); medemp1[1]=sumc(seqa(1,1,m).*mfre1[.,1]); @ Variable2 @ mm=dati[.,2]; dd=saved(mm,"comodo2","m2"); call openvars("comodo2"); {bmat1,bmat2} = CUB("comodo2","m2",fatt1,fatt1); output off; cmlset; pigre0[2]=1/(1+exp(-bmat1[1])); pigre1[2]=1/(1+exp(-bmat1[1]-bmat1[2])); CSI0[2]=1/(1+exp(-bmat1[3])); CSI1[2]=1/(1+exp(-bmat1[3]-bmat1[4])); mteo0[.,2]=PROBCUB00(m,pigre0[2],csi0[2]); mteo1[.,2]=PROBCUB00(m,pigre1[2],csi1[2]); mfre0[.,2]=counts(dati0[.,2],seqa(1,1,m))/rows(dati0); mfre1[.,2]=counts(dati1[.,2],seqa(1,1,m))/rows(dati1); med0[2]=sumc(seqa(1,1,m).*mteo0[.,2]); med1[2]=sumc(seqa(1,1,m).*mteo1[.,2]); medemp0[2]=sumc(seqa(1,1,m).*mfre0[.,2]); medemp1[2]=sumc(seqa(1,1,m).*mfre1[.,2]); freqdop0=freqdoppia(dati0[.,1],dati0[.,2]); freqdop1=freqdoppia(dati1[.,1],dati1[.,2]); F0=cumsumc(mteo0[.,2]); G0=cumsumc(mteo0[.,1]); F1=cumsumc(mteo1[.,2]); G1=cumsumc(mteo1[.,1]); matfre=(freqdop0*rows(dati0))~(freqdop1*rows(dati1)); _cml_Bounds=(0~100)|(0~100); _cml_CovPar=1; _cml_Diagnostic=1; {pa,b1,b2,b3,b4}=cml(matfre,0,&likpsicov,1.2|0.2); cmlset; {HHMAT0,PP0}=Hpro(F0,G0,exp(pa[1])); {HHMAT1,PP1}=Hpro(F1,G1,exp(pa[1]+pa[2])); output on; "________________________________________"; "________________________________________"; "average of CUB distribution fitted to colour"; "cova=0" med0[1]; "cova=1" med1[1]; "sample mean"; "cova=0" medemp0[1]; "cova=1" medemp1[1]; "________________________________________"; "CSI cova=0" csi0[1]; "CSI cova=1" csi1[1]; "________________________________________"; " pigre gamma0 gamma1"; mat1[1:3]'; "standard error"; sqrt(diag(mat2))'; "________________________________________"; "________________________________________"; "average of CUB distribution fitted to flavour "; "cova=0" med0[2]; "cova=1" med1[2]; "sample mean"; "cova=0" medemp0[2]; "cova=1" medemp1[2]; "________________________________________"; "PI cova=0" pigre0[2]; "PI cova=1" pigre1[2]; "________________________________________"; "CSI cova=0" csi0[2]; "CSI cova=1" csi1[2]; "________________________________________"; "beta0 beta1 gamma0 gamma1"; bmat1[1:4]'; "standard error"; sqrt(diag(bmat2))'; "________________________________________"; "log-likel=" likePAE; "ETA" pa'; "psi"; "cova=0" exp(pa[1]); "cova=1" exp(pa[1]+pa[2]); "________________________________________"; "probab (X>5, y>5)"; "Lower knowl.= " sumc(sumc(pp0[6:7,6:7])); "High knowl.= " sumc(sumc(pp1[6:7,6:7])); "________________________________________"; "________________________________________"; output off; @::::::JACKKNIFE::::::::::::@ __m=7; m=7; nboot=rows(dati); mteo0s=zeros(m,2); mteo1s=zeros(m,2); mfre0s=zeros(m,2); mfre1s=zeros(m,2); csi0s=zeros(2,1); pigres=zeros(2,1); csi1s=zeros(2,1); psis=zeros(2,1); ncofin=rows(pigres|csi0s|csi1s)+5; finalej=zeros(nboot,ncofin); codice=zeros(nboot,1); riga=1; do while riga<=nboot; if riga==nboot-1; datis=dati[2:rows(dati),.]; fatt1S=datis[.,3]; elseif riga==nboot; datis=dati[1:rows(dati)-1,.]; fatt1S=datis[.,3]; else; datis=dati[1:riga riga+2:rows(dati),.]; fatt1S=datis[.,3]; endif; dati0s=SELIF(datis,fatt1S.==0); dati1s=SELIF(datis,fatt1S.==1); @Variable1 @ mm=datis[.,1]; dd=saved(mm,"comodo1","m1"); call openvars("comodo1"); cmlset; {smat1, smat2}=CUB("comodo1","m1",0,fatt1S); pigres[1]=smat1[1]; CSI0s[1]=1/(1+exp(-smat1[2])); CSI1s[1]=1/(1+exp(-smat1[2]-smat1[3])); mteo0s[.,1]=PROBCUB00(m,pigres[1],csi0s[1]); mteo1s[.,1]=PROBCUB00(m,pigres[1],csi1s[1]); mfre0s[.,1]=counts(dati0s[.,1],seqa(1,1,m))/rows(dati0s); mfre1s[.,1]=counts(dati1s[.,1],seqa(1,1,m))/rows(dati1s); @Variable2@ cmlset; __m=7; m=7; mm=datis[.,2]; dd=saved(mm,"comodo2","m2"); call openvars("comodo2"); cmlset; {sbmat1,sbmat2} = CUB("comodo2","m2",0,fatt1S); pigres[2]=sbmat1[1]; CSI0s[2]=1/(1+exp(-sbmat1[2])); CSI1s[2]=1/(1+exp(-sbmat1[2]-sbmat1[3])); mteo0s[.,2]=PROBCUB00(m,pigres[2],csi0[2]); mteo1s[.,2]=PROBCUB00(m,pigres[2],csi1[2]); mfre0s[.,2]=counts(dati0s[.,2],seqa(1,1,m))/rows(dati0s); mfre1s[.,2]=counts(dati1s[.,2],seqa(1,1,m))/rows(dati1s); freqdop0s=freqdoppia(dati0s[.,1],dati0s[.,2]); freqdop1s=freqdoppia(dati1s[.,1],dati1s[.,2]); F0=cumsumc(mteo0s[.,2]); G0=cumsumc(mteo0s[.,1]); F1=cumsumc(mteo1s[.,2]); G1=cumsumc(mteo1s[.,1]); matfres=(freqdop0s*rows(dati0s))~(freqdop1s*rows(dati1s)); cmlset; _cml_Bounds=(0~100)|(0~100); _cml_CovPar=1; _cml_Diagnostic=1; {pas,b1s,b2s,b3s,b4s}=cml(matfres,0,&likpsicov,pa); codice[riga]=b4s; cmlset; finalej[RIGA,.]=pigres'~csi0s'~csi1s'~pas'~likepae~exp(pas[1])~exp(sumc(pas)); riga=riga+1; endo; output on; "__________________________________________________"; "JACKNIFE: ESTIMATED standard error of psi"; "COVA=0 COVA=1"; (stdc(finalej[.,10 11])*sqrt(((nboot-1)^2)/nboot))'; output off; end;