sbgcop

sbgcop: Semiparametric Bayesian Gaussian copula estimation This package estimates parameters of a G

.packageName <- "sbgcop"
"ldmvnorm" <-
function(Y,S) {     # log density of a matrix with mvn rows
 n<-dim(Y)[1]
 p<-dim(Y)[2]
 -.5*n*log(det(S)) -.5*n*p*log(2*pi)-.5*sum( diag( solve(S)%*%t(Y)%*%Y)) }

"plot.psgc" <-
function(x,...){ 

Y<-x$Y.pmean
p<-dim(Y)[2]
vnames<-colnames(Y)
nc<-max( nchar(vnames))


par(mar=c(.6,2,.5,1),mgp=c(.6,.5,0),oma=c(.75*nc,2.85,1,1),xpd=FALSE,cex=.8)
par(mfcol=c(p,3))
###
for(j in 1:p){
y<-Y[,j]
if(length(unique(y))<=10 ){
barplot(table(y)/sum(!is.na(y)),xlab="",ylab=vnames[j],xaxt="n",yaxt="n")
                           }
if(length(unique(y))>10 ){
hist(y,prob=TRUE,main="",xlab="",ylab=vnames[j],,xaxt="n",yaxt="n",col="gray")
                          }
              }
###
plotci.sA(x$C.psamp,ylabs=rep("",p),mgp=c(.6,.5,0))
plotci.sA(sR.sC(x$C.psamp),ylabs=rep("",p),mgp=c(.6,.5,0))
###
                         }

"plotci.sA" <-
function(sA,ylabs=colnames(sA[,,1]),mgp=c(1.75,.75,0)) {
qA<-qM.sM(sA)
p<-dim(qA)[1]
tmp<-c(qA)
tmp<-tmp[tmp!=1]
par(mgp=mgp)

for(j in 1:p) {
plot(0,0,type="n",ylim=range(c(tmp),na.rm=TRUE),xlim=c(1,p),
     ylab=ylabs[j],xaxt="n", xlab="")
points( (1:p)[-j], qA[j,-j,2],pch=16,cex=.6 )
segments( x0=(1:p)[-j], y0=qA[j,-j,1], x1=(1:p)[-j], y1=qA[j,-j,3] )
abline(h=0,col="gray")
abline(v=j,col="gray")
               }
axis(side=1,at=1:p, labels=colnames(qA[,,1]),las=2)
                                              }

"print.sum.psgc" <-
function(x,...) {
cat("\n #### MCMC details         ####","\n")
cat("\n number of saved samples:",x$nsamp,"\n")
cat("\n average effective sample size:",mean(x$ESS),"\n")
cat("\n effective sample sizes","\n" )
print(round(x$ESS,2))

cat("\n #### Parameter estimation ####","\n")
cat(" Posterior quantiles of correlation coefficients:","\n")
print(round(x$QC,2))

cat("\n Posterior quantiles of regression coefficients:","\n")
print(round(x$QR,2))


                                }

"qM.sM" <-
function(sM,quantiles=c(0.025,.5,.975)) {

p1<-dim(sM)[1]
p2<-dim(sM)[2]
s <-dim(sM)[3]

qM<-array(dim=c(p1,p2,length(quantiles)))
dimnames(qM)<-list(dimnames(sM)[[1]],dimnames(sM)[[2]], 
              paste(quantiles*100,rep("% quantile",length(quantiles)),sep=""))

for(l in 1:length(quantiles)) {
qM[,,l]<-apply(sM,c(1,2),quantile,prob=quantiles[l],na.rm=TRUE)
                               }
qM
}

"rwish" <-
function(S0,nu){       # sample from a Wishart distribution
sS0<-chol(S0)
Z<-matrix(rnorm(nu*dim(S0)[1]),nu,dim(S0)[1])%*%sS0
t(Z)%*%Z                     }

"sR.sC" <-
function(sC) {
p<-dim(sC)[1]
s<-dim(sC)[3]
sR<-array(dim=c(p,p,s) )
dimnames(sR)<-dimnames(sC)

for(l in 1:s) {

C<-sC[,,l]
R<-C*NA
for(j in 1:p){
R[j,-j]<- C[j,-j]%*%solve(C[-j,-j])
             }
sR[,,l]<-R     
              }
sR                    }

"sbgcop.mcmc" <-
function(Y,S0=diag(dim(Y)[2]),n0=dim(Y)[2]+2, 
                      nsamp=100,odens=max(1,round(nsamp/1000)),seed=1, 
                      verb=TRUE){

########## check input
ok_S0<-all(eigen(S0)$val>0) & dim(S0)[1]==dim(Y)[2] & dim(S0)[2]==dim(Y)[2]
ok_n0<-(n0>=0)

if(!ok_S0) { cat("Error: S0 must be a positive definite p x p matrix \n") }
if(!ok_n0) { cat("Error: n0 must be positive \n") }

if(ok_S0 & ok_n0) {

########## data
n<-dim(Y)[1]
p<-dim(Y)[2]
##########

########## starting values
set.seed(seed)
R<-NULL
for(j in 1:p) { R<-cbind(R, match(Y[,j],sort(unique(Y[,j])))) }
Rlevels<-apply(R,2,max,na.rm=TRUE)
Ranks<- apply(Y,2,rank,ties.method="max",na.last="keep")
N<-apply(!is.na(Ranks),2,sum)
U<- t( t(Ranks)/(N+1))
Z<-qnorm(U)
Zfill<-matrix(rnorm(n*p),n,p)
Z[is.na(Y)]<-Zfill[is.na(Y) ]
S<-cov(Z)
##########

########## things to keep track of
Y.pmean<-matrix(0,nrow=n,ncol=p)
LPC<-NULL
C.psamp<-array(dim=c(p,p,floor(nsamp/odens)))
dimnames(C.psamp)<-list(colnames(Y),colnames(Y),1:floor(nsamp/odens))
##########

########## start MCMC
for(ns in 1:nsamp) {

#### update Z
for(j in sample(1:p)) {
Sjc<- S[j,-j]%*%solve(S[-j,-j])
sdj<- sqrt( S[j,j] -S[j,-j]%*%solve(S[-j,-j])%*%S[-j,j]  )
muj<- Z[,-j]%*%t(Sjc)

for(r in 1:Rlevels[j]){

ir<- (1:n)[R[,j]==r & !is.na(R[,j])]
lb<-suppressWarnings(max( Z[ R[,j]==r-1,j],na.rm=TRUE))
ub<-suppressWarnings(min( Z[ R[,j]==r+1,j],na.rm=TRUE))
Z[ir,j]<-qnorm(runif(length(ir),
               pnorm(lb,muj[ir],sdj),pnorm(ub,muj[ir],sdj)),muj[ir],sdj)
                       }
ir<-(1:n)[is.na(R[,j])]
Z[ir,j]<-rnorm(length(ir),muj[ir],sdj)
                          }
#### 

#### update S
S<-solve(rwish(solve(S0*n0+t(Z)%*%Z),n0+n))
####

#### save results
if(ns %% odens ==0) {
C<-S/(sqrt(diag(S))%*%t(sqrt(diag(S))))

lpc<-ldmvnorm(Z%*%diag(1/sqrt(diag(S))),C)
LPC<-c(LPC,lpc)
C.psamp[,,ns/odens]<-C


Y.imp<-Y
for(j in 1:p) {
Y.imp[is.na(Y[,j]),j]<-quantile(Y[,j],
                       pnorm(Z[is.na(Y[,j]),j],0,sqrt(S[j,j])),na.rm=TRUE,type=1)
                }
Y.pmean<-  ( (ns/odens-1)/(ns/odens) )*Y.pmean+  (1/(ns/odens) )*Y.imp
                  }
if(verb==TRUE & (ns %%(odens*10))==0){
        cat(round(100*ns/nsamp),"percent done ",date(),"\n")
                                     }
####

         }
##########

G.ps<-list(C.psamp=C.psamp,Y.pmean=Y.pmean,LPC=LPC)
class(G.ps)<-"psgc"
return(G.ps)
              }
}

"summary.psgc" <-
function(object,...){
qC<-qM.sM(object$C.psamp)
qR<-qM.sM(sR.sC(object$C.psamp))

p<-dim(qC)[1]
nsamp<-dim(object$C.psamp)[3]
vn<-colnames(qC[,,1])
ACF<-QC<-QR<-NULL

rnamesC<-rnamesR<-NULL

for(l1 in 1:p)  {
for(l2 in (1:p)[-l1]) {

QR<-rbind(QR, c(qR[l1,l2,1:3]) )
rnamesR<-c(rnamesR,paste(vn[l1],vn[l2],sep="~")  )

if(l1<l2) {
QC<-rbind(QC, c(qC[l1,l2,1:3]) )
ACF<-rbind(ACF,acf(object$C.psamp[l1,l2,],lag.max=round(nsamp/20),plot=FALSE)[[1]][-1] )
rnamesC<-c(rnamesC, paste(vn[l1],vn[l2],sep="*")  ) 
            }
 }}

rownames(QC)<-rownames(ACF)<-rnamesC
rownames(QR)<-rnamesR

Kappa<-1+2*apply(ACF,1,sum)

ESS<-nsamp/Kappa

ACR.lag1<-matrix(0,p,p)
ACR.lag1[upper.tri(ACR.lag1)]<-ACF[1:choose(p,2),1]
ACR.lag1<-ACR.lag1+t(ACR.lag1)
diag(ACR.lag1)<-NA

nsamp<-dim(object$C.psamp)[3]

SUM<-list(QC=QC,QR=QR,nsamp=nsamp,ACR.lag1=ACR.lag1,ESS=ESS)
structure(SUM,class="sum.psgc")

}