Nothing
# includes functions: compAllocVtoCompAlloc, ls.mvrm, label.count, plot.generic, plot.mvrm
compAllocVtoCompAlloc <- function(G,p,compAllocV){
compAlloc<-array(-19,p*(p-1)/2)
move <- 0
for (k in 0:(G-1)){
temp <- 0
for (i in 1:(p-1)){
for (j in (i+1):p){
if (compAllocV[i] == compAllocV[j] && compAllocV[i] == k){
compAlloc[(i-1)*p+j-(i)*(i+1)/2] <- move
temp <- temp+1
}
}
}
if (temp > 0) move<-move+1
}
move2<-0
for (k in 0:(G-2)){
for (l in (k+1):(G-1)){
temp <- 0
for (i in 1:(p-1)){
for (j in (i+1):p){
if (((compAllocV[i] == k && compAllocV[j] == l) || (compAllocV[i] == l && compAllocV[j] == k))){
compAlloc[(i-1)*p+j-(i)*(i+1)/2] <- move+move2
temp<-temp+1
}
}
}
if (temp > 0) move2<-move2+1
}
}
return(compAlloc)
}
ls.mvrm <- function(x){
file1 <- paste(x$DIR,"BNSP.nu.ls.txt",sep="")
file2 <- paste(x$DIR,"BNSP.eta.ls.txt",sep="")
file3 <- paste(x$DIR,"BNSP.nmembers.ls.txt",sep="")
file4 <- paste(x$DIR,"BNSP.clusters.txt",sep="")
nu<-mvrm2mcmc(x,"nu")
eta<-mvrm2mcmc(x,"eta")
nmembers<-mvrm2mcmc(x,"nmembers")
nparams<-x$LGc+(x$LGc+1)+1 #gamma + eta + nmembers
nu <- nu[,c(mapply(seq,seq(1,x$LGc),length.out=x$H,by=x$LGc))]
eta <- eta[,c(mapply(seq,seq(1,x$LGc+1),length.out=x$H,by=x$LGc+1))]
mat<-cbind(nu,eta,nmembers)
mcmc.pars <- array(data = mat, dim = c(x$nSamples, x$H, nparams))
compAlloc<-mvrm2mcmc(x,"compAlloc")
compAllocV<-mvrm2mcmc(x,"compAllocV")
ifelse(x$G == 1, z<-compAlloc+1, z<-compAllocV+1)
probs<-mvrm2mcmc(x,"probs")
ifelse(x$G ==1, DIM<-x$H, DIM<-x$G)
pr <- array(data = probs, dim = c(x$d, x$nSamples, DIM))
#pr <- array(data = probs, dim = c(x$p, x$nSamples, x$H))
p <- aperm(pr, c(2,1,3))
if (x$G == 1){
ls <- label.switching(method = "STEPHENS", p = p, mcmc = mcmc.pars, z = z, K = x$H)
mcmc.pars.per<-permute.mcmc(mcmc.pars, ls$permutations$STEPHENS)$output
}
if (x$G > 1){
ls <- label.switching(method = "STEPHENS", p = p, z = z, K = x$G)
#Permute component allocations for the variables
compAllocVPer<-compAllocV
for (sw in 1:x$nSamples)
for (k in 0:(x$G-1))
compAllocVPer[sw,compAllocV[sw,]==ls$permutations[[1]][sw,k+1]-1]<-k
#Permute component allocations for the correlations
compAllocPer<-matrix(0,ncol=x$d,nrow=x$nSamples)
for (sw in 1:x$nSamples){
if (sum(compAllocVPer[sw,]==compAllocV[sw,])==x$p) compAllocPer[sw,]<-compAlloc[sw,]
else compAllocPer[sw,]<-compAllocVtoCompAlloc(x$G,x$p,compAllocVPer[sw,])
}
#Find clustering
tab <- table(apply(compAllocPer+1, 1, paste, collapse=", "))
ls$clusters <- as.numeric(unlist(strsplit(names(which.max(tab)), ', ')))
#From permutation matrix of variables to permuation matrix of the correlations
fpm<-matrix(-99,nrow=x$nSamples,ncol=x$H)
for (sw in 1:x$nSamples){
b<-compAlloc[sw,]
a<-compAllocPer[sw,]
fpv<-array(-99,x$H)
if ((max(a,b)+2)<=x$H) fpv[(max(a,b)+2):x$H]<-c((max(a,b)+2):x$H)
for (k in 0:max(a,b))
fpv[k+1]<-a[min(which(b==k))]+1
fpm[sw,] <- fpv
}
#With fpm, permute parameters
mcmc.pars.per<-permute.mcmc(mcmc.pars, fpm)$output
}
nu<-NULL
if (x$LGc > 0){
nu<-mcmc.pars.per[,,1:x$LGc]
nu <-matrix(nu,nrow=x$nSamples)
nu <- nu[,c(mapply(seq,seq(1,x$H),length.out=x$LGc,by=x$H))]
}
eta<-mcmc.pars.per[,,(1+x$LGc):(nparams-1)]
eta <-matrix(eta,nrow=x$nSamples)
eta <- eta[,c(mapply(seq,seq(1,x$H),length.out=(x$LGc+1),by=x$H))]
nmembers<-mcmc.pars.per[,,nparams]
write.table(nu,file=file1,row.names = FALSE,col.names = FALSE)
write.table(eta,file=file2,row.names = FALSE,col.names = FALSE)
write.table(nmembers,file=file3,row.names = FALSE,col.names = FALSE)
write.table(ls$clusters,file=file4,row.names = FALSE,col.names = FALSE, quote=FALSE)
}
label.count <- function(term,labels,counts,model){
if (is.character(term)){
int.label1<-int.label2<-grep(term,labels,fixed=TRUE)
k<-1
while (!length(int.label1)==1 && (nchar(term)>k)){
int.label1<-grep(substr(term,1,nchar(term)-k),labels,fixed=TRUE)
k<-k+1
}
k<-1
while (!length(int.label2)==1 && (nchar(term)>k)){
int.label2<-grep(substr(term,k,nchar(term)),labels,fixed=TRUE)
k<-k+1
}
if ((!length(int.label1) == 1) || (!length(int.label2) == 1) || !(int.label1==int.label2))
stop(cat(cat("`term` in", model, "model should be an integer between 1 and", length(labels),"or one of: "), cat(labels,sep=", ",fill=TRUE)),call. = FALSE)
int.label<-int.label1
}
if (is.numeric(term)) int.label<-term
if ((!length(int.label) == 1) || (int.label > length(labels)))
stop(cat(cat("`term` in", model, "model should be an integer between 1 and", length(labels),"or one of: "), cat(labels,sep=", ",fill=TRUE)),call. = FALSE)
count<-counts[int.label]
return(c(int.label,count))
}
plotGeneric <- function(mvrmObj, MEAN, STDEV, CORM, CORS, DEP, SCALE,
response, intercept, grid,
centre, quantiles, static, contour, centreEffects, plotOptions,
int.label, count, vars, label, is.D, formula.term, which.Spec, assign,
knots, storeMeanVector, storeIndicator, data, plotEmptyCluster){
p<-mvrmObj$p
period0<-0
for (i in 1:length(formula.term)){
if (!grepl("knots",formula.term[i]) && which.Spec[i] > 0){
formula.term[i]<-substr(formula.term[i],1,nchar(formula.term[i])-1)
formula.term[i]<-paste(formula.term[i],", knots=knots[[",which.Spec[i],"]])",sep="")
}
if (grepl("period = 0",formula.term[i])){
period0<-1
calc.p <- mean(mvrm2mcmc(mvrmObj,"period"))
formula.term[i]<-substr(formula.term[i],1,nchar(formula.term[i])-1)
formula.term[i]<-sub(", period = 0", "", formula.term[i])
formula.term[i]<-paste(formula.term[i],", period = ", calc.p,")",sep="")
}
}
V<-NULL
Dynamic<-NULL
for (k in 1:length(int.label)) V<-c(V,which(assign %in% int.label[k]))
if (STDEV) V <- V - 1
if ((intercept || sum(is.D)) && MEAN) V<-c(1,V)
indeces<-V
if (MEAN && !(1 %in% indeces)) indeces<-c(1,indeces)
if (STDEV) indeces<-c(1,indeces+1)
if (count==1){
if (!is.D){
max.shift <- 0
if (mvrmObj$nBreaks > 0) max.shift <- max(mvrm2mcmc(mvrmObj,"breaks"))
min1<-min(with(data,eval(as.name(vars[1]))))
max1<-max(with(data,eval(as.name(vars[1])))) + max.shift
newR1<-seq(min1,max1,length.out=grid)
newData<-data.frame(newR1)
colnames(newData)<-vars
if (sum(!is.na(formula.term))){
dm<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=knots,
meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)
Ds<-dm$X
Dynamic<-dm$Dynamic
amplitude<-dm$DSP[1]
harmonics<-dm$DSP[3]
startSin<-dm$DSP[2]
ifelse(intercept, startSin <- startSin, startSin <- startSin - 1)
}
if (!sum(!is.na(formula.term))) Ds<-matrix(1,1,1)
if (((!1%in%V)||STDEV) && sum(!is.na(formula.term))) Ds<-Ds[,-1]
}
if (is.D){
newData<-data.frame(sort(unique(with(data,eval(as.name(vars[1]))))))
colnames(newData)<-vars
Ds<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=NULL,
meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)$X
if (STDEV) Ds<-Ds[,-1]
}
}
if (count==2){
if (sum(is.D)){
dv<-which(is.D==1)
cv<-which(is.D==0)
min1<-min(with(data,eval(as.name(vars[cv]))))
max1<-max(with(data,eval(as.name(vars[cv]))))
newR1<-seq(min1,max1,length.out=grid)
newR2<-unique(with(data,eval(as.name(vars[dv]))))
newData<-expand.grid(newR1,newR2)
colnames(newData)<-vars[c(cv,dv)]
#which.Spec <- which.Spec[which.Spec > 0]
#whichKnots <- which.Spec
#Dstar<-knots[[whichKnots]]
#Ds<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=Dstar,
# meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)$X
Ds<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=knots,
meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)$X
if (STDEV) Ds<-Ds[,-1]
}
if (sum(is.D)==0){
min1<-min(with(data,eval(as.name(vars[1]))))
min2<-min(with(data,eval(as.name(vars[2]))))
max1<-max(with(data,eval(as.name(vars[1]))))
max2<-max(with(data,eval(as.name(vars[2]))))
newR1<-seq(min1,max1,length.out=grid)
newR2<-seq(min2,max2,length.out=grid)
newData<-expand.grid(newR1,newR2)
colnames(newData)<-vars
#whichKnots <- which.Spec
#Dstar<-knots[[whichKnots]]
#Ds<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=Dstar,
# meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)$X
Ds<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=knots,
meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)$X
if (STDEV) Ds<-Ds[,-1]
}
}
if (MEAN){
dim2<-1
if (CORM) dim2<-mvrmObj$H
#print("jbf")
#print(dim(Ds))
fit<-array(0,dim=c(mvrmObj$nSamples,dim2,NROW(Ds)))
#fit<-matrix(0,nrow=mvrmObj$nSamples,ncol=NROW(Ds))
#meanReg<-array(0,dim=c(x$nSamples,x$H,x$LUT))
if (CORM==0 && DEP==0){
etaFN <- file.path(paste(mvrmObj$DIR,"BNSP.beta.txt",sep=""))
how.many <- p*(mvrmObj$LG+1)
multiplier <- mvrmObj$LG+1
if (!is.null(Dynamic)){
harmonicsFN <- file.path(paste(mvrmObj$DIR,"BNSP.Hbeta.txt",sep=""))
how.many2 <- 2*harmonics
}
if (period0==1){
periodFN <- file.path(paste(mvrmObj$DIR,"BNSP.period.txt",sep=""))
periodFile<-file(periodFN,open="r")
per.val <- calc.p
}
}else if (CORM==1){
if (mvrmObj$H==1) etaFN <- file.path(paste(mvrmObj$DIR,"BNSP.eta.txt",sep=""))
if (mvrmObj$H > 1){
etaFN <- file.path(paste(mvrmObj$DIR,"BNSP.eta.ls.txt",sep=""))
clusFN <- paste(mvrmObj$DIR,"BNSP.clusters.txt",sep="")
if (!file.exists(etaFN) || !file.exists(clusFN))
quiet(ls.mvrm(mvrmObj))
clus<-array(unlist(read.table(clusFN)))
tabs<-table(clus)
labs<-array(0,mvrmObj$H)
labs[as.numeric(names(tabs))]<-tabs
V <- rep(V,mvrmObj$H) + rep(seq(0,max(V)*(mvrmObj$H-1),by=max(V)),each=length(V))
}
how.many <- (mvrmObj$LGc+1)*mvrmObj$H
multiplier<-0
}else{
etaFN <- file.path(paste(mvrmObj$DIR,"BNSP.psi.txt",sep=""))
how.many <- p*p*mvrmObj$LK
multiplier <- mvrmObj$LK
}
eFile<-file(etaFN,open="r")
if (!is.null(Dynamic))
eFile2<-file(harmonicsFN,open="r")
for (i in 1:mvrmObj$nSamples){
eta<-scan(eFile,what=numeric(),n=how.many,quiet=TRUE)
if (!is.null(Dynamic)){
etaHar<-scan(eFile2,what=numeric(),n=how.many2,quiet=TRUE)
ifelse(period0==1,
per.val<-scan(periodFile,what=numeric(),n=1,quiet=TRUE),
per.val<-mvrmObj$period)
sinWphs <- 0 * newR1
for (l in 1:harmonics)
sinWphs <- sinWphs +
etaHar[2*l-1] * sin(2 * l * pi * newR1 / per.val) +
etaHar[2*l] * cos(2 * l * pi * newR1 / per.val)
Ds[,c((startSin+1):(startSin+amplitude+1))] <- sweep(Dynamic, MARGIN=1, sinWphs, `*`)
}
if (period0==1 && is.null(Dynamic)){
per.val.r <- per.val
per.val<-scan(periodFile,what=numeric(),n=1,quiet=TRUE)
#for (k in 1:length(formula.term)){
# formula.term[k]<-substr(formula.term[k],1,nchar(formula.term[k])-1)
# formula.term[k]<-sub(paste(", period = ", per.val.r,sep=""), "", formula.term[k])
# formula.term[k]<-paste(formula.term[k],", period = ", per.val,")",sep="")
#}
#dm<-DM(formula=reformulate(formula.term),data=newData,n=NROW(newData),knots=knots,
# meanVector=storeMeanVector[indeces],indicator=storeIndicator[indeces],centre=TRUE)
#Ds<-dm$X
for (l in 1:harmonics)
Ds[,startSin+2*l-1] <- sin(2 * l * pi * newR1 / per.val)
Ds[,startSin+2*l] <- cos(2 * l * pi * newR1 / per.val)
}
fit[i,,] <- matrix(c(eta[V+(response-1)*multiplier]),byrow=TRUE,ncol=dim(as.matrix(Ds))[2]) %*% t(as.matrix(Ds))
if (CORM && mvrmObj$FT==1) fit[i,,] <- tanh(fit[i,,])
if (centreEffects) fit[i,,]<-fit[i,,]-mean(fit[i,,])
}
close(eFile)
if (!is.null(Dynamic)) close(eFile2)
if (period0==1) close(periodFile)
}
if (STDEV){
fit<-array(0,dim=c(mvrmObj$nSamples,1,NROW(Ds)))
if (CORS==0 && SCALE==0){
alphaFN <- file.path(paste(mvrmObj$DIR,"BNSP.alpha.txt",sep=""))
sigma2FN <- file.path(paste(mvrmObj$DIR,"BNSP.sigma2.txt",sep=""))
how.many1 <- p*mvrmObj$LD
how.many2 <- p
}
if (CORS==1){
alphaFN <- file.path(paste(mvrmObj$DIR,"BNSP.omega.txt",sep=""))
sigma2FN <- file.path(paste(mvrmObj$DIR,"BNSP.sigma2R.txt",sep=""))
how.many1 <- mvrmObj$LDc
how.many2 <- 1
}
if (SCALE==1){
alphaFN <- file.path(paste(mvrmObj$DIR,"BNSP.psi.txt",sep=""))
sigma2FN <- file.path(paste(mvrmObj$DIR,"BNSP.phi2.txt",sep=""))
how.many1 <- p*mvrmObj$LC
how.many2 <- p
}
aFile<-file(alphaFN,open="r")
s2File<-file(sigma2FN,open="r")
for (i in 1:mvrmObj$nSamples){
alpha<-scan(aFile,what=numeric(),n=how.many1,quiet=TRUE)
#ifelse(intercept || !is.na(formula.term),
# s2<-scan(s2File,what=numeric(),n=how.many2,quiet=TRUE)[response], s2<-1)
#if (!is.na(formula.term) && intercept)
# fit[i,,]<-s2*exp(as.matrix(Ds)%*%matrix(c(alpha[V+(response-1)*(how.many1/p)])))
#if (!is.na(formula.term) && !intercept)
# fit[i,,]<-exp(as.matrix(Ds)%*%matrix(c(alpha[V+(response-1)*(how.many1/p)])))
#if (is.na(formula.term)) fit[i,,]<-s2
ifelse(intercept, s2<-scan(s2File,what=numeric(),n=how.many2,quiet=TRUE)[response], s2<-1)
fit[i,,]<-s2*exp(as.matrix(Ds)%*%matrix(c(alpha[V+(response-1)*(how.many1/p)])))
if (!SCALE) fit[i,,]<-sqrt(fit[i,,])
if (centreEffects) fit[i,,]<-fit[i,,]/mean(fit[i,,])
}
close(aFile)
close(s2File)
}
if (count==1 && !is.D){
newX<-unlist(newData)
centreM<-drop(apply(fit,c(2,3),centre))
if (sum(!is.na(formula.term)))
dataM<-data.frame(rep(newX,dim(fit)[2]),c(t(centreM)))
if (!sum(!is.na(formula.term)))
dataM<-data.frame(expand.grid(newX,centreM))
nms <- c(vars,"centreM")
if (!is.null(quantiles)){
nms<-c(nms,"QLM","QUM")
QM<-apply(fit,c(2,3),quantile,probs=quantiles,na.rm=TRUE)
if (sum(!is.na(formula.term)))
dataM<-data.frame(dataM,c(t(QM[1,,])),c(t(QM[2,,])))
if (!sum(!is.na(formula.term)))
dataM<-data.frame(dataM,rep(QM[1,,],each=length(newX)),rep(QM[2,,],each=length(newX)))
}
colnames(dataM) <- nms
if (CORM==1 && mvrmObj$H > 1){
dataM<-data.frame(dataM,group=factor(rep(1:mvrmObj$H,each=grid)),size=rep(labs,each=grid))
dataM<-dataM[order(-dataM$size, dataM[,1]),]
labs2<-labs
if (!plotEmptyCluster) {
dataM<-dataM[dataM$size > 0,]
labs2<-labs[labs>0]
}
}
#dataRug<-data.frame(b=with(data,eval(as.name(vars))))
plotElM<-c(geom_line(aes_string(x=as.name(vars),y=centreM),col=4,alpha=0.5),
#geom_rug(mapping=aes(x=dataRug$b),data=dataRug,alpha=0.3),
list(ylab(label)))
if (!is.null(quantiles))
plotElM<-c(geom_ribbon(data=dataM,aes_string(x=vars, ymin="QLM", ymax="QUM"),alpha=0.2),plotElM)
if (CORM==1 && mvrmObj$H > 1){
plotElM<-c(geom_line(aes_string(x=as.name(vars),y="centreM",group="group",linetype="group"),col=4,alpha=0.5),
#geom_rug(mapping=aes(x=dataRug$b),data=dataRug,alpha=0.3),
list(ylab(label)))
if (!is.null(quantiles))
plotElM<-c(geom_ribbon(data=dataM,aes_string(x=vars,ymin="QLM",ymax="QUM",group="group"),alpha=0.2),plotElM)
}
ggM<-ggplot(data=dataM)
plotM<-ggM + plotElM + plotOptions
if (CORM==1 && mvrmObj$H > 1){
cns<-noquote(unlist(lapply(seq(1,p-1),function (k) paste(rep(k,p-k),seq(k+1,p),sep=""))))
cgsa<-noquote(paste(lapply(1:length(labs2), function(k)cns[clus==k])))
plotM <- ggM + plotElM + guides(fill=FALSE) +
scale_linetype_manual(name = "correlation groups: ",
labels = cgsa,
values=seq(1:length(labs2))) +
theme(legend.position = "bottom") +
plotOptions
}
}
if (count==1 && is.D){
lvs<-levels(with(data,eval(as.name(vars[1]))))
if (is.null(lvs)) lvs<-sort(unique(with(data,eval(as.name(vars[1])))))
#sort(unique(with(data,eval(as.name(vars[1])))))
df<-data.frame(x=rep(lvs,each=mvrmObj$nSamples),y=c(fit))
plotElM<-c(geom_boxplot(),list(xlab(label),ylab("")))
#ggM<-ggplot(data=df,aes(x=factor(df$x),y=df$y))
ggM<-ggplot(df,aes(factor(df[,1]),df[,2]))
# ggM<-ggplot(data=dataM,aes_string(x=vars[1],y=vars[2],z=centreM))
plotM<-ggM + plotElM + plotOptions
}
if (count==2 && sum(is.D)==1){
centreM<-drop(apply(fit,c(2,3),centre))
disc.var<-vars[which(is.D==1)]
cont.var<-vars[which(is.D==0)]
dataM<-data.frame(newData,centreM)
nms <- c(colnames(newData),"centreM")
if (!is.null(quantiles)){
nms<-c(nms,"QLM","QUM")
QM<-apply(fit,c(2,3),quantile,probs=quantiles,na.rm=TRUE)
dataM<-data.frame(dataM,c(t(QM[1,,])),c(t(QM[2,,])))
}
colnames(dataM) <- nms
#DG<-data.frame(b=with(data,eval(as.name(cont.var))),c=with(data,eval(as.name(disc.var))))
#plotElM<-c(geom_line(aes_string(x=as.name(vars),y="centreM",group="group",linetype="group"),col=4,alpha=0.5),
# geom_rug(mapping=aes(x=dataRug$b),data=dataRug,alpha=0.3),
# list(ylab(label)))
#if (!is.null(quantiles))
# plotElM<-c(geom_ribbon(data=dataM,aes_string(x=vars,ymin="QLM",ymax="QUM",group="group"),alpha=0.2),plotElM
plotElM<-c(geom_line(aes_string(x=cont.var,y=centreM,
group=disc.var,linetype=disc.var,col=disc.var),alpha=0.5),
#geom_rug(mapping=aes(x=DG$b,group=c,linetype=c),data=DG,alpha=0.3),
list(ylab(label[which(is.D==0)])))
if (!is.null(quantiles))
plotElM<-c(geom_ribbon(data=dataM,aes_string(x=cont.var,ymin="QLM",ymax="QUM",
group=disc.var,col=disc.var),alpha=0.2),plotElM)
ggM<-ggplot(data=dataM)
plotM<-ggM + plotElM + plotOptions
}
if (count==2 && sum(is.D)==0){
centreM<-drop(apply(fit,c(2,3),centre))
if (contour==1){
dataM<-data.frame(newData,centreM)
nms <- c(colnames(newData),"centreM")
colnames(dataM) <- nms
ggM<-ggplot(data=dataM,aes_string(x=vars[1],y=vars[2],z=centreM))
level<-1
plotM<-ggM + geom_contour(data=dataM,aes(colour = stat(level))) + plotOptions
}
if (static && contour==0){
defaultList<-list(x=as.numeric(newR1),y=as.numeric(newR2),z=matrix(centreM,length(newR1),length(newR2)),colvar=matrix(centreM,length(newR1),length(newR2)))
along="xy";
space=0.6;
if (MEAN) optionalList<-list(xlab=vars[1],ylab=vars[2],zlab=label,along=along,space=space,add=FALSE,bty="g",main=paste("mean of", mvrmObj$varsY[[response]]))
if (STDEV) optionalList<-list(xlab=vars[1],ylab=vars[2],zlab=label,along=along,space=space,add=FALSE,bty="g",main=paste("stdev of", mvrmObj$varsY[[response]]))
allOptions<-c(defaultList,plotOptions,optionalList)
do.call("ribbon3D",allOptions[!duplicated(names(allOptions))])
}
if (!static && contour==0){
a<-as.matrix(cbind(newData,centreM))
if (is.null(plotOptions$col)) col=rainbow(16,2/3)
else col<-plotOptions$col
plotCentreM <- centreM
if (min(centreM)<=0) plotCentreM <- centreM + abs(min(centreM)) + 1
ra<-ceiling(length(col)*plotCentreM/max(plotCentreM))
defaultList<-list(x=a,col=col[ra])
if (MEAN) optionalList<-list(size=0.4,bg=1,axisLabels=c(vars[1],label,vars[2]),main=paste("mean of", mvrmObj$varsY[[response]]))
if (STDEV) optionalList<-list(size=0.4,bg=1,axisLabels=c(vars[1],label,vars[2]),main=paste("stdev of", mvrmObj$varsY[[response]]))
allOptions<-c(defaultList,plotOptions,optionalList)
do.call("scatterplot3js",allOptions[!duplicated(names(allOptions))])
}
}
if (contour==1)
return(plotM)
}
plot.mvrm <- function(x, model, term, response, response2,
intercept = TRUE, grid = 30, centre = mean,
quantiles = c(0.1, 0.9), contour = TRUE, static = TRUE,
centreEffects = FALSE, plotOptions = list(), nrow, ask = FALSE,
plotEmptyCluster = FALSE, combine = FALSE, ...)
{
oldpar <- NULL
on.exit(par(oldpar))
#
if (!is.function(centre)) stop("centre must be a function, usually mean or median")
if (length(quantiles)==1) quantiles <- c(quantiles,1-quantiles)
if (length(quantiles) > 2) stop("up to two quantiles")
if (!is.null(quantiles)) quantiles <- unique(sort(quantiles))
grid<-round(grid)
if (missing(response)) response <- c(1:x$p)
if (missing(response2)) response2 <- c(1:x$p)
if (max(response) > x$p) stop("argument response exceeds the number of responses")
#
if (missing(model)) model<-"all"
MEAN <- 0; STDEV <- 0; CORM <- 0; CORS <- 0; DEP <- 0; SCALE <- 0
if ((model=="all" || model=="mean") && (x$NG > 0)) MEAN <- 1
if ((model=="all" || model=="stdev") && (x$ND > 0)) STDEV <- 1
if ((model=="all" || model=="corm") && x$p > 1 && x$LUT > 1) CORM <- 1
if ((model=="all" || model=="cors") && x$p > 1 && x$LUT > 1) CORS <- 1
if ((model=="all" || model=="dep") && (x$NK > 0)) DEP <- 1
if ((model=="all" || model=="scale") && (x$NC > 0)) SCALE <- 1
#
if (MEAN==0 && STDEV==0 && CORM==0 && CORS==0 && DEP==0 && SCALE==0) stop("no terms to plot");
termM <- NULL; termSD <- NULL; termMC<-NULL; termSDC<-NULL; termD<-NULL; termSC<-NULL
countM <- NULL; countSD <- NULL; countMC <- NULL; countSDC <- NULL; countD <- NULL; countSC <- NULL
int.labelM <- NULL; int.labelSD <- NULL; int.labelD <- NULL; int.labelSC <- NULL
if (missing(term)){
if (MEAN) {termM<-1:x$NG; countM<-x$countX; int.labelM <- 1:x$NG}
if (STDEV) {termSD<-1:x$ND; countSD<-x$countZ; int.labelSD <- 1:x$ND}
if (DEP) {termD<-1:(x$NK-1); countD<-x$countC[termD]; int.labelD <- termD}
if (SCALE) {termSC<-1:x$NC; countSC<-x$countW; int.labelSC <- 1:x$NC}
}else{
if (MEAN){
termM<-term
for (i in termM){
#print(i)
lcX <- label.count(i,x$labelsX,x$countX,"mean")
int.labelM[i]<-lcX[1]
countM[i]<-lcX[2]
}
}
if (STDEV){
termSD<-term
for (i in termSD){
lcZ <- label.count(i,x$labelsZ,x$countZ,"stdev")
int.labelSD[i]<-lcZ[1]
countSD[i]<-lcZ[2]
}
}
if (DEP){
termD<-term
for (i in termD){
lcC <- label.count(i,x$labelsC,x$countC,"autoregressive")
int.labelD[i]<-lcC[1]
countD[i]<-lcC[2]
}
}
if (SCALE){
termSC<-term
for (i in termSC){
lcW <- label.count(i,x$labelsW,x$countW,"scale")
int.labelSC[i]<-lcW[1]
countSC[i]<-lcW[2]
}
}
}
if (CORM) {termMC<-1; countMC<-x$countXc[termMC]}
if (CORS) {termSDC<-1; countSDC<-x$countZc[termSDC]}
#
if ((contour == 0) &&
(sum(c(countM,countSD,countMC,countSDC,countD,countSC)>1) > 0) &&
(length(c(countM,countSD,countMC,countSDC,countD,countSC)) > 1)){
warning("for 3-d plots, specify one model, one term and one response at a time")
contour <- 1
}
#
if (combine && length(term) == 1) combine <- FALSE
if (combine) termM <- termSD <- DEP <- SCALE <- 1
#
my_plots <- list()
count <- 1
for (r in response){
if (MEAN)
for (i in termM){
#lcX <- label.count(i,x$labelsX,x$countX,"mean")
#int.labelX <- lcX[1]
#countX <- lcX[2]
if (!combine){
int.labelX <- int.labelM[i]
countX <- countM[i]
varsArg <- x$varsX[[int.labelX]]
is.DArg <- x$is.Dx[[int.labelX]]
which.SpecArg <- x$which.SpecX[[int.labelX]]
}else{
int.labelX <- term
varsArg <- unique(unlist(x$varsX[int.labelX]))
countX <- length(varsArg)
is.DArg <- unlist(x$is.Dx[int.labelX])
if (!sum(is.DArg)) is.DArg <- 0
which.SpecArg <- unlist(x$which.SpecX[int.labelX])
if ((length(unique(varsArg)) > 1) && (sum(!is.DArg) > 1)) stop("can't combine more than 2 continuous terms; consider interaction terms")
}
if (countX==1) contour <- 1
plotOptions2 <- list(ggtitle(paste("mean of", x$varsY[r])),plotOptions)
my_plots[[count]] <- plotGeneric(mvrmObj=x, MEAN=1, STDEV=0, CORM=0, CORS=0, DEP=0, SCALE=0,
response=r, intercept=intercept, grid=grid,
centre=centre, quantiles=quantiles, static=static,
contour=contour,
centreEffects=centreEffects, plotOptions=plotOptions2,
int.label=int.labelX, count=countX, vars=varsArg,
label=x$labelsX[int.labelX], is.D=is.DArg,
formula.term=x$formula.termsX[int.labelX],
which.Spec=which.SpecArg, assign=x$assignX,
knots=x$Xknots, storeMeanVector=x$storeMeanVectorX,
storeIndicator=x$storeIndicatorX, data=x$data,
plotEmptyCluster=plotEmptyCluster)
if ((ask==TRUE) && (length(my_plots)>0)) print(my_plots[[count]])
if (count==1)
oldpar <- c(oldpar, par(ask=ask))
count <- count + 1
}
if (STDEV)
for (i in termSD){
#lcZ <- label.count(i,x$labelsZ,x$countZ,"stdev")
#int.labelZ <- lcZ[1]
#countZ <- lcZ[2]
if (!combine){
int.labelZ <- int.labelSD[i]
countZ <- countSD[i]
varsArg <- x$varsZ[[int.labelZ]]
is.DArg <- x$is.Dz[[int.labelZ]]
which.SpecArg <- x$which.SpecZ[[int.labelZ]]
}else{
int.labelZ <- term
countZ <- 2
varsArg <- unlist(x$varsZ[int.labelZ])
is.DArg <- unlist(x$is.Dz[int.labelZ])
which.SpecArg <- unlist(x$which.SpecZ[int.labelZ])
}
if (combine && (!sum(is.DArg)==1)) stop("can only combine a continuous with a discrete term")
if (countZ==1) contour <- 1
plotOptions2 <- list(ggtitle(paste("stdev of", x$varsY[[r]])),plotOptions)
if (x$LUT > 1) plotOptions2 <- list(ggtitle(paste("innov stdev of", x$varsY[[r]])),plotOptions)
my_plots[[count]] <- plotGeneric(mvrmObj=x, MEAN=0, STDEV=1, CORM=0, CORS=0, DEP=0, SCALE=0,
response=r, intercept=intercept, grid=grid,
centre=centre, quantiles=quantiles, static=static,
contour=contour,
centreEffects=centreEffects, plotOptions=plotOptions2,
int.label=int.labelZ, count=countZ, vars=varsArg,
label=x$labelsZ[int.labelZ], is.D=is.DArg,
formula.term=x$formula.termsZ[int.labelZ],
which.Spec=which.SpecArg, assign=x$assignZ,
knots=x$Zknots, storeMeanVector=x$storeMeanVectorZ,
storeIndicator=x$storeIndicatorZ, data=x$data,
plotEmptyCluster=plotEmptyCluster)
if ((ask==TRUE) && (length(my_plots)>0)) print(my_plots[[count]])
if (count==1)
oldpar <- c(oldpar, par(ask=ask))
count <- count + 1
}
if (SCALE)
for (i in termSC){
#lcW <- label.count(i,x$labelsW,x$countW,"scale")
#int.labelW <- lcW[1]
#countW <- lcW[2]
if (!combine){
int.labelW <- int.labelSC[i]
countW <- countSC[i]
varsArg <- x$varsW[[int.labelW]]
is.DArg <- x$is.Dw[[int.labelW]]
which.SpecArg <- x$which.SpecW[[int.labelW]]
}else{
int.labelW <- term
countW <- 2
varsArg <- unlist(x$varsW[int.labelW])
is.DArg <- unlist(x$is.Dw[int.labelW])
which.SpecArg <- unlist(x$which.SpecW[int.labelW])
}
if (combine && (!sum(is.DArg)==1)) stop("can only combine a continuous with a discrete term")
if (countW==1) contour <- 1
plotOptions2 <- list(ggtitle(paste("scale of", x$varsY[[r]])),plotOptions)
my_plots[[count]] <- plotGeneric(mvrmObj=x, MEAN=0, STDEV=1, CORM=0, CORS=0, DEP=0, SCALE = 1,
response=r, intercept=intercept, grid=grid,
centre=centre, quantiles=quantiles, static=static,
contour=contour,
centreEffects=centreEffects, plotOptions=plotOptions2,
int.label=int.labelW, count=countW, vars=varsArg,
label=x$labelsW[int.labelW], is.D=is.DArg,
formula.term=x$formula.termsW[int.labelW],
which.Spec=which.SpecArg, assign=x$assignW,
knots=x$Wknots, storeMeanVector=x$storeMeanVectorW,
storeIndicator=x$storeIndicatorW, data=x$data,
plotEmptyCluster=plotEmptyCluster)
if ((ask==TRUE) && (length(my_plots)>0)) print(my_plots[[count]])
if (count==1)
oldpar <- c(oldpar, par(ask=ask))
count <- count + 1
}
}
if (DEP)
for (r1 in response){
for (r2 in response2){
Pair<-(r1-1)*x$p + r2
for (i in termD){
#lcC <- label.count(i,x$labelsC,x$countC,"autoregressive")
#int.labelC <- lcC[1]
#countC <- lcC[2]
if (!combine){
int.labelC <- int.labelD[i]
countC <- countD[i]
varsArg <- x$varsC[[int.labelC]]
is.DArg <- x$is.Dc[[int.labelC]]
which.SpecArg <- x$which.SpecC[[int.labelC]]
}else{
int.labelC <- term
countC <- 2
varsArg <- unlist(x$varsC[int.labelC])
is.DArg <- unlist(x$is.Dc[int.labelC])
which.SpecArg <- unlist(x$which.SpecC[int.labelC])
}
if (combine && (!sum(is.DArg)==1)) stop("can only combine a continuous with a discrete term")
if (countC==1) contour <- 1
plotOptions2 <- list(ggtitle(paste("autoreg", x$varsY[r1], x$varsY[r2])), plotOptions)
my_plots[[count]] <- plotGeneric(mvrmObj=x, MEAN=1, STDEV=0, CORM=0, CORS=0, DEP=1, SCALE=0,
response=Pair, intercept=intercept, grid=grid,
centre=centre, quantiles=quantiles, static=static,
contour=contour,
centreEffects=centreEffects, plotOptions=plotOptions2,
int.label=int.labelC, count=countC, vars=varsArg,
label=x$labelsC[int.labelC], is.D=is.DArg,
formula.term=x$formula.termsC[int.labelC],
which.Spec=which.SpecArg, assign=x$assignC,
knots=x$Cknots, storeMeanVector=x$storeMeanVectorC,
storeIndicator=x$storeIndicatorC, data=x$lag,
plotEmptyCluster=plotEmptyCluster)
if ((ask==TRUE) && (length(my_plots)>0)) print(my_plots[[count]])
if (count==1)
oldpar <- c(oldpar, par(ask=ask))
count <- count + 1
}
}
}
if (CORM){
for (i in termMC){
contour <- 1
plotOptions2 <- list(ggtitle(paste("mean cor")),plotOptions)
which.Spec <- 0
if (length(x$which.SpecXc)) which.Spec <- x$which.SpecXc[[1]]
is.D <- 0
if (length(x$is.Dxc)) is.D <- x$is.Dxc[[1]]
vars <- x$varTime
if (length(x$varsXc)) vars <- x$varsXc[[1]]
my_plots[[count]] <- plotGeneric(mvrmObj=x, MEAN=1, STDEV=0, CORM=1, CORS=0, DEP=0, SCALE=0,
response=1, intercept=intercept,
grid=grid, centre=centre, quantiles=quantiles, static=static, contour=contour,
centreEffects=centreEffects, plotOptions=plotOptions2,
int.label=1, count=1, vars=vars,
label=x$labelsXc[1], is.D=is.D,
formula.term=x$formula.termsXc[1],
which.Spec=which.Spec, assign=x$assignXc,
knots=x$Xcknots, storeMeanVector=x$storeMeanVectorXc,
storeIndicator=x$storeIndicatorXc, data=x$data,
plotEmptyCluster=plotEmptyCluster)
if ((ask==TRUE) && (length(my_plots)>0)) print(my_plots[[count]])
if (count==1)
oldpar <- c(oldpar, par(ask=ask))
count <- count + 1
}
}
if (CORS){
for (i in termSDC){
contour <- 1
plotOptions2 <- list(ggtitle(paste("stdev cor")),plotOptions)
which.Spec <- 0
if (length(x$which.SpecZc)) which.Spec <- x$which.SpecZc[[1]]
is.D <- 0
if (length(x$is.Dzc)) is.D <- x$is.Dzc[[1]]
vars <- x$varTime
if (length(x$varsZc)) vars <- x$varsXc[[1]]
my_plots[[count]] <- plotGeneric(mvrmObj=x, MEAN=0, STDEV=1, CORM=0, CORS=1, DEP=0, SCALE=0,
response=1, intercept=intercept,
grid=grid, centre=centre, quantiles=quantiles, static=static, contour=contour,
centreEffects=centreEffects, plotOptions=plotOptions2,
int.label=1, count=1, vars=vars,
label=x$labelsZc[1], is.D=is.D,
formula.term=x$formula.termsZc[1],
which.Spec=which.Spec, assign=x$assignZc,
knots=x$Zcknots, storeMeanVector=x$storeMeanVectorZc,
storeIndicator=x$storeIndicatorZc, data=x$data,
plotEmptyCluster=plotEmptyCluster)
if ((ask==TRUE) && (length(my_plots)>0)) print(my_plots[[count]])
if (count==1)
oldpar <- c(oldpar, par(ask=ask))
count <- count + 1
}
}
if (missing(nrow)){
if (CORM==0 && CORS==0 && DEP==0){nrow <- min(x$p,length(my_plots))}else{nrow <- ceiling(sqrt(count-1))}
}
if ((ask==FALSE) && (length(my_plots)>0)) quiet(print(grid.arrange(grobs = my_plots, nrow = nrow)))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.