## Functions used in the intermiss.R library(statmod) ######################################################################## ##### Functions to fit MIssing Data Models, Run Genetic Algorithms etc. ###################################### ########################################################################### fun.y1<- function(y1, muimis, covyimis,gimis, gpimis,eta10imis,eta20imis,eta11imis, gamma1) { n=length(y1) res= numeric(n) py1= numeric(n) mu1=muimis[1] sd1= sqrt(covyimis[1,1]) var21= covyimis[2,2]- covyimis[1,2]^2/covyimis[1,1] mu21= muimis[2] + covyimis[1,2]/covyimis[1,1]*(y1-mu1) py1= fun.probmiss(gimis[1],gpimis[1],eta10imis[1],eta20imis[1],eta11imis[1], y1,gamma1) for (i in 1:length(y1)) { res[i]=integrate(fun.nllnigimiss, lower=-Inf, upper=Inf, nimis=1, muimis=mu21[i], covyimis=var21, gimis=gimis[2], gpimis=gpimis[2], eta10=eta10imis[2],eta20=eta20imis[2],eta11=eta11imis[2],gamma1=gamma1)$value } res= res * py1 * dnorm(y1, mean=mu1, sd=sd1) res } fun.csmat<-function(sigma,rho,nobs,transform,case=1) { ## calculate the the related quantity for compound symmetry vcov matrix. ## input ## sigma -- the standard error of the error term ## rho --- the correlation coefficient ## nobs --- the number of observations, nobs>=1 ## transform --- logical number indicator if the parameters are transformed ## output ## case=1 --- nobs x nobs CS matrix ## case=2 --- nobs x nobs derivative matrix wrt sigma ## case=3 --- nobs x nobs derivative matrix wrt rho if (transform==FALSE) { if (case==1) res<-sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs)) ## derivative wrt to sigma else if (case==2) res<- 2*sigma*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs)) ## derivative wrt to rho else if (case==3) res<- sigma^2*(matrix(1,nrow=nobs,ncol=nobs)-diag(1,nrow=nobs)) } else { sigma<-fun.trans2pos(sigma) rho <- fun.trans2rho(rho) if (case==1) res<- sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs)) else if (case==2) res<- 2*sigma^2*(rho*matrix(1,nrow=nobs,ncol=nobs)+diag(1-rho,nrow=nobs)) else if (case==3) res<- sigma^2*(matrix(1,nrow=nobs,ncol=nobs)-diag(1,nrow=nobs))* (1-rho^2)/2 } res } fun.table<-function(data, time, mtype=1) { ## function to generate data for ISNI calculation ## input ## data: panel dataset in the long format, include "id", "time", "y" ## row for per subject at ALL timepoints, no matter missing or not. ## time: all the planned time visits. ## mtype: 1 -- three missingness type ## 2 -- only two missingness type. ## output ## a dataset from first time to the dropout time or to the last time point ## more variables added: ## yp -- previous observed outcome ## g -- missingness indicator, 0-observed, 1-intermittent missing, 2-dropout ## gp -- missingness indicator in the previous visit. ## first sort the data by id and time data = data[order(data$id, data$time),] uid = unique(data$id) lastT=numeric(length(uid)) maxT=max(time) minT=min(time) yp=g=gp=res=NULL if (mtype ==1 ) { for (i in 1:length(uid)) { print(uid[i]) datai = data[data$id==uid[i],,drop=F] ypi =gi=gpi=numeric(dim(datai)[1]) lastT[i] = max(datai$time[!is.na(datai$y)]) if (lastT[i] != maxT) lastT[i] = min(time[time>lastT[i]]) for (j in 1:dim(datai)[1]) { ## assign values to yp, g, gp if (datai$time[j] == minT) { ypi[j]=0;gi[j]= -1; gpi[j]=-1 } else { prevT = max(time[timelastT[i]]) for (j in 1:dim(datai)[1]) { ## assign values to yp, g, gp if (datai$time[j] == minT) { ypi[j]=0;gi[j]= -1; gpi[j]=-1 } else { prevT = max(time[time2 & data$gp==0 & data$gp2==0,] if (nrow(data00)>0) { if (any(data00$g==1) & any(data00$g==2)){ mdm00=multinom(gmodelyr57, data=data00, trace=FALSE, Hess=TRUE,maxit=500) data00=cbind(data00, fitprobs=mdm00$fit[,1])} else if (any(data00$g==1) | any(data00$g==2)) { mdm00=multinom(gmodelyr57, data=data00, trace=FALSE, Hess=TRUE,maxit=500) data00=cbind(data00, fitprobs=1-mdm00$fit[,1])} else {mdm00=NULL; data00= cbind(data00, fitprobs=1)} } data01= data[data$time >2 & data$gp==0 & data$gp2==1,] if (nrow(data01)>0) { if (any(data01$g==1) & any(data01$g==2)){ mdm01=multinom(gmodelyr57, data=data01, trace=FALSE, Hess=TRUE,maxit=500) data01=cbind(data01, fitprobs=mdm01$fit[,1])} else if (any(data01$g==1) | any(data01$g==2)) { mdm01=multinom(gmodelyr57, data=data01, trace=FALSE, Hess=TRUE,maxit=500) data01=cbind(data01, fitprobs=1-mdm01$fit[,1])} else {mdm01=NULL; data01= cbind(data01, fitprobs=1)} } data10= data[data$time >2 & data$gp==1 & data$gp2==0,] if (nrow(data10)>0) { if (any(data10$g==1) & any(data10$g==2)){ mdm10=multinom(gmodelyr57, data=data10, trace=FALSE, Hess=TRUE,maxit=500) data10=cbind(data10, fitprobs=mdm10$fit[,1])} else if (any(data10$g==1) | any(data10$g==2)) { mdm10=multinom(gmodelyr57, data=data10, trace=FALSE, Hess=TRUE,maxit=500) data10=cbind(data10, fitprobs=1-mdm10$fit[,1])} else {mdm10=NULL; data10= cbind(data10, fitprobs=1)} } data11= data[data$time >2 & data$gp==1 & data$gp2==1,] if (nrow(data11)>0) { if (any(data11$g==1) & any(data11$g==2)){ mdm11=multinom(gmodelyr57, data=data11, trace=FALSE, Hess=TRUE,maxit=500) data11=cbind(data11, fitprobs=mdm11$fit[,1])} else if (any(data11$g==1) | any(data11$g==2)) { mdm11=multinom(gmodelyr57, data=data11, trace=FALSE, Hess=TRUE,maxit=500) data11=cbind(data11, fitprobs=1-mdm11$fit[,1])} else {mdm11=NULL; data11= cbind(data11, fitprobs=1)} } data.mdm = rbind(data.vm1,data.yr2,data00, data01, data10, data11) data.mdm = data.mdm[order(data.mdm$id, data.mdm$time),] res=list(data=data.mdm, mdm.yr2=mdm.yr2, mdm00=mdm00, mdm01=mdm01, mdm10=mdm10, mdm11=mdm11) return(res) } else stop("No lag more than 2") } fun.gfit2=function(data,gmodel, lag=2){ ## missing data model fitting when all the missingness are deemed to be potential transitory. ## Therefore all the missingness are intermittent. ## first create lag-2 missingness indicator data=data[,!names(data)=="fitprobs"] if (any(names(data)=="fitprobs")) stop("There should be no fitprobs in the data") gmodel0=gmodel$gmodel0 ## model at intermittent visit with gp=0 gmodel1=gmodel$gmodel1 ## model at intermittent visit with gp=1 firstT= min(data$time) lastT = max(data$time) ## now fit models if (lag==1) { ## first, baseline visit0 data.vm1 = cbind(data[data$time==firstT,], fitprobs=1) ## then for year 2, 5, and 7 gmodel= g~black+female+blackfemale+year5+year7+blackt5+blackt7 data0=data[data$time>0 & data$gp==0,] if (nrow(data0)>0) { if (any(data0$g==1)) {mdm0=glm(gmodel, data=data0, family=binomial) data0=cbind(data0, fitprobs=1-mdm0$fit)} else {mdm0=NULL; data0= cbind(data0, fitprobs=1)} } data1=data[data$time>0 & data$gp==1,] if (nrow(data1)>0) { if (any(data1$g==1)) {mdm1=glm(gmodel, data=data1, family=binomial) data1=cbind(data1, fitprobs=1-mdm1$fit)} else {mdm1=NULL; data1= cbind(data1, fitprobs=1)} } data.mdm = rbind(data.vm1,data0, data1) data.mdm = data.mdm[order(data.mdm$id, data.mdm$time),] res=list(data=data.mdm, mdm0=mdm0,mdm1=mdm1) return(res) } if (lag==2) { ## first, baseline visit 0 data.vm1 = cbind(data[data$time==firstT,], fitprobs=1) ## then, for visit at year 2 gmodelyr2= g~black+female+blackfemale+yp data.yr2= data[data$time==2,] if (any(data.yr2$g==1)) {mdm.yr2= glm(gmodelyr2, data=data.yr2, family=binomial) data.yr2=cbind(data.yr2, fitprobs=1-mdm.yr2$fit)} else {mdm.yr2=NULL; data.yr2=cbind(data.year2, fitprobs=1)} ## then, for visit at year 5, and 7 mdm00=mdm01=mdm10=mdm11=NULL gmodelyr57= g~black+female +blackfemale+year7+blackt7+yp data00= data[data$time >2 & data$gp==0 & data$gp2==0,] if (nrow(data00)>0) { if (any(data00$g==1)) {mdm00=glm(gmodelyr57, data=data00, family=binomial) data00=cbind(data00, fitprobs=1-mdm00$fit)} else {mdm00=NULL; data00= cbind(data00, fitprobs=1)} } data01= data[data$time >2 & data$gp==0 & data$gp2==1,] if (nrow(data01)>0) { if (any(data01$g==1)) {mdm01=glm(gmodelyr57, data=data01, family=binomial) data01=cbind(data01, fitprobs=1-mdm01$fit)} else {mdm01=NULL; data01= cbind(data01, fitprobs=1)} } data10= data[data$time >2 & data$gp==1 & data$gp2==0,] if (nrow(data10)>0) { if (any(data10$g==1)) {mdm10=glm(gmodelyr57, data=data10, family=binomial) data10=cbind(data10, fitprobs=1-mdm10$fit)} else {mdm10=NULL; data10= cbind(data10, fitprobs=1)} } data11= data[data$time >2 & data$gp==1 & data$gp2==1,] if (nrow(data11)>0) { if (any(data11$g==1)) {mdm11=glm(gmodelyr57, data=data11, family=binomial) data11=cbind(data11, fitprobs=1-mdm11$fit)} else {mdm11=NULL; data11= cbind(data11, fitprobs=1)} } data.mdm = rbind(data.vm1,data.yr2,data00, data01, data10, data11) data.mdm = data.mdm[order(data.mdm$id, data.mdm$time),] res=list(data=data.mdm, mdm.yr2=mdm.yr2, mdm00=mdm00, mdm01=mdm01, mdm10=mdm10, mdm11=mdm11) return(res) } } fun.dermean4gen = function(gendata,gvec, pglw.ig, fix.form, rdm.form,r1grp){ ## generate the derivative of mean for genetic algorithm. ## g in data cannot be 2. all the missingness should be coded as 1. gendata=gendata[order(gendata$ambig, gendata$id, decreasing=T),] gendata$gvec=0 gendata$gvec[1:(4*length(gvec))]=rep(gvec,each=4) gendata$g=ifelse(gendata$ambig==0, gendata$g, ifelse(gendata$gvec==1, gendata$ig, gendata$dg)) gendata$gp=ifelse(gendata$ambig==0, gendata$gp, ifelse(gendata$gvec==1, gendata$igp, gendata$dgp)) gendata$gp2=ifelse(gendata$ambig==0, gendata$gp2, ifelse(gendata$gvec==1, gendata$igp2, gendata$dgp2)) if(any(gendata$g==2)) stop("no dropout (g=2) allowed") gen.dmean=fun.glmmall(pglw.ig$par,data=gendata,fix.form=yfix.form,rdm.form=yrdm.form1,r1grp=r1grp, transform=FALSE, case=4) return(gen.dmean) } fun.gen<-function(gvec,gendata, pglw.nabla11, para, gen.dmean, yfix.form, yrdm.form1, r1group, slope=1) { ## first generate dataset with g imputed by gvec ##gendata=pglw gendata=gendata[order(gendata$ambig, gendata$id, decreasing=T),] gendata$gvec=0 gendata$gvec[1:(4*length(gvec))]=rep(gvec,each=4) gendata$g=ifelse(gendata$ambig==0, gendata$g, ifelse(gendata$gvec==1, gendata$ig, gendata$dg)) gendata$gp=ifelse(gendata$ambig==0, gendata$gp, ifelse(gendata$gvec==1, gendata$igp, gendata$dgp)) gendata$gp2=ifelse(gendata$ambig==0, gendata$gp2, ifelse(gendata$gvec==1, gendata$igp2, gendata$dgp2)) ## for (i in 1:length(uid)) { ## print(i) ## datai=pglw[pglw$id==uid[i],,drop=F] ## if (datai$ambig[1]==1 ) { ## if (gvec[j]==1) { ## datai$g=datai$ig; datai$gp=datai$igp; datai$gp2=datai$igp2 ## } else ## if (gvec[j]==2){ ## datai$g=datai$dg; datai$gp= datai$dgp; datai$gp2=datai$dgp2 ## } else ## stop("no such gvec number") ## j=j+1 ## idedata[idedata$id==uid[i],]=datai ## } ## ## } ## if any of G is dropout if (any(gendata$g==2)) gen.gres=fun.gfit(data=gendata, gmodel=gmodel, lag=2) else ## for only two missingness type gen.gres= fun.gfit2(data=gendata, gmodel=gmodel,lag=2) gendata=gen.gres$data gendata=gendata[order(gendata$ambig, gendata$id, decreasing=T),] ## calculate the ISNI temp = matrix(0, nrow=length(gen.dmean), ncol=length(para)) uid=unique(gendata$id) for (i in 1:length(uid)) { if (!is.null(gen.dmean[[i]])){ ##datai=gendata[gendata$id==uid[i],,drop=F] ##gfitimis <- datai$fitprobs[is.na(datai$y)] gfitimis= ##gendata$fitprobs[((i-1)*4+1):(4*i)][is.na(gendata$y[((i-1)*4+1):(4*i)])] gendata$fitprobs[gendata$id==uid[i] & is.na(gendata$y)] temp[i,]=gen.dmean[[i]][,1:length(gfitimis)]%*% as.matrix(gfitimis) } } ur1= unique(r1group) ur1= ur1[order(ur1)] gen.nabla12=matrix(0, nrow=length(para), ncol=length(ur1)) for (j in 1:length(ur1)) { tempj=temp[r1group==ur1[j],,drop=F] gen.nabla12[,j] = apply(tempj,2,sum) } gen.isni= pglw.nabla11%*%gen.nabla12 gen.isni= abs(gen.isni[,1])+abs(gen.isni[,2])+ abs(gen.isni[,3]) + abs(gen.isni[,4]) print(gen.isni) return(gen.isni[slope]) } ######################################################################## ##### ISNI for GLMM ###################################### ########################################################################### fun.glmmall <- function (par, data,fix.form, rdm.form, transform=TRUE,r1grp=NULL, subdiv=100, nqpt=10, case=1) { ## This function is used to fit the ignorable GLMM, and obtain \nabla_{11} and \nabla_{12} ## Input ## par : 1 x (nbeta+nD) vector for the parameter$\beta$and$D$to be maximized ## data : data.frameof the entire data (for y, yp, g and gp), including observed data and ## observations at intermittent and dropout time. ## it must have the following components ## id, time, y, yp, g, gp ## Note: the value of gp should never be 2 ## The data should be otained from the fun.table ## columns for fixed effect covariates in y, randome effects for y ## predictors for g. ## fix.form: formular of the fixed effect model ## rdm.form: formular of the random effect model, y~1 is the random intercept model. ## case : list of items to be calcuated ## 1: A scalar negative marginal loglikelihood is to be calculated ## 2: A vector of the first derivative of negative marginal loglikelihood w.r.t ##$\beta$and$D$. ## 3: A vector of the first derivative of the conditional mean of missing data ## given the observed data w.r.t$\beta$and$D$. ## 4: return a list of der.meanymi ## r1grp --- names of the variables for the grouping indicator for \gamma_1 vector ## nqpt --- number of quadrature point if (any(data$gp==2)) stop("The value of gp should not be 2") ##nmiss = length(unique(data$id[data$g>0])) ##dmeanym=vector("list",length=nmiss) fix.temp<-model.frame(fix.form,data=data,na.action=NULL) x<- as.matrix(model.matrix(fix.form, data=fix.temp)) y<- as.vector(model.extract(fix.temp, "response")) rdm.temp<-model.frame(rdm.form,data=data, na.action=NULL) z<- as.matrix(model.matrix(rdm.form, data=rdm.temp)) nbeta<-dim(x)[2] beta<- par[1:nbeta] D = par[(nbeta+1):length(par)] if (length(D) != dim(z)[2]) stop("The number of random effects is not compatible") if (length(par) != (dim(x)[2]+dim(z)[2]) ) stop("Length of Parameter is incorrect") uid = unique(data$id) if (case==1 ) temp = matrix(0, nrow=length(uid), ncol=1) if (case==2 | case==3) temp = matrix(0, nrow=length(uid), ncol=length(par)) if (case==4) temp=vector("list", length=length(uid)) if (!is.null(r1grp)) r1group=numeric(length(uid)) for (i in 1:length(uid)) { xi= as.matrix(x[data$id==uid[i],,drop=F]) zi= as.matrix(z[data$id==uid[i],,drop=F]) yi= y[data$id==uid[i]] datai = list(yi=yi,xi=xi,zi=zi) ## print(xi) if (!is.null(r1grp)) r1group[i]=data[data$id==uid[i],,drop=F][1,r1grp] if (case <3) { ##temp[i,] = fun.glmmsubiadapt(datai,beta=beta,D=D, transform=transform, case=case,subdiv=subdiv) temp[i,] = fun.glmmsubi(datai,beta=beta,D=D, transform=transform,nqpt=nqpt, case=case) }else if (case==3) { gfiti= data$fitprobs[data$id==uid[i]] if (any(is.na(yi))) { temp[i,] = fun.glmmsubi(datai, beta=beta,D=D, transform=transform, gfiti=gfiti, nqpt=nqpt,case=case) } } else ## for case 4 { if(any(is.na(yi))) { temp[[i]] = fun.glmmsubi(datai, beta=beta,D=D, transform=transform, nqpt=nqpt,case=case) } }##end of ifelse loop } ## end of for loop if (case==1) res= sum(temp) if (case==2) res= apply(temp, 2,sum) if (case==3) { if (is.null(r1grp)) res= apply(temp, 2, sum) else { ur1= unique(r1group) ur1= ur1[order(ur1)] res=matrix(0, nrow=length(par), ncol=length(ur1)) for (j in 1:length(ur1)) { tempj=temp[r1group==ur1[j],,drop=F] ## print(tempj) res[,j] = apply(tempj,2,sum) } } } if (case==4) res=list(dmean=temp, r1group=r1group) res } fun.glmmsubi <- function ( datai,beta, D,transform, gfiti=NULL, nqpt, case=1) { ## This function is to calculate the quantities in the subject level ## input ## datai -- a list of data for the ith subject ## yi -- ni x 1 vector of the response for the ith subject ## xi -- ni x nbeta matrix of the covariates for the fixed effect ## zi -- ni x nD matrix of the covariates for the random effect ## beta -- 1 x nbeta vector of parameters for the fixed effect # D -- parameters in the variance-covariance matrix for the random effects ## transform -- wether or not the parameter in the D is transformed. ## case -- 1, calculate the negative loglikelihood for subject i ## 2, calculate the negative score value for subject i ## 3, calcuate the \nabla_12 for subject i. ## 4, calcuate the der.meanymi ## output ## res -- if case=1 , res=negative loglikelihood ## if case=2, res=negative score value ## if case=3, res=\nabla_12 for subject i. ## if case=4, res=der.meanymi ## obtain marginal likelihood for subject i. if (case==1) mar.lik.subi<-fun.ghqd(fun.subiq,datai,nqpt,beta,D,transform=transform,case=1) else if ( case ==2 ) { ## obtain derivative of marginal loglikelihood for subject i mar.lik.subi<-fun.ghqd(fun.subiq,datai,nqpt,beta,D,transform=transform,case=1) der.mlogl.subi <-fun.ghqd(fun.subiq,datai,nqpt,beta,D,transform=transform,case=2)/mar.lik.subi } else if ( (case==3 |case==4) & any(is.na(datai$yi)) ){ gfitimis <- gfiti[is.na(datai$yi)] mar.lik.subi<-fun.ghqd(fun.subiq,datai,nqpt,beta,D,transform=transform,case=1) ##print(mar.lik.subi) der.mlik.subi<-fun.ghqd(fun.subiq,datai,nqpt,beta,D,transform=transform,case=2) der.meanymi = (fun.ghqd(fun.dmean1,datai,nqpt, beta,D,transform=transform)/mar.lik.subi +fun.ghqd(fun.dmean2,datai,nqpt,beta,D,transform=transform)/mar.lik.subi -as.matrix(der.mlik.subi) %*% t(as.matrix(fun.ghqd(fun.dmean3,datai,nqpt,beta,D,transform=transform))) /(mar.lik.subi^2)) ##print(gfitimis) ## print(der.meanymi) if (case==3) nabla12i= der.meanymi %*% as.matrix(gfitimis) } if (case==1) res<- -log(mar.lik.subi) ## A scalar negative marginal loglikelihood if (case==2) {res<- -der.mlogl.subi} ## derivative of the negative marginal logliklihood if (case==3) res<- nabla12i ##der.meanymi if (case==4) res= der.meanymi res } fun.dmean1<-function(qpt,datai,beta,D,transform) fun.subiq(qpt,datai,beta,D,transform=transform,case=3)$part1 fun.dmean2<-function(qpt,datai,beta,D,transform) fun.subiq(qpt,datai,beta,D,transform=transform,case=3)$part2 fun.dmean3<-function(qpt,datai,beta,D,transform) fun.subiq(qpt,datai,beta,D,transform=transform,case=3)$part3 fun.subiq<- function (qpt,datai,beta,D,transform,case=1) { ## Fun.subiq is to calculate the following items (listed in case) for subject i at quadrature point qpt. ## qpt : A scalar quadrature point ## yi -- ni x 1 vector of the response for the ith subject ## xi -- ni x nbeta matrix of the covariates for the fixed effect ## zi -- ni x nD matrix of the covariates for the random effect ## beta -- 1 x nbeta vector of parameters for the fixed effect ## D -- parameters in the variance-covariance matrix for the random effects: standard deviation ## transform -- wether or not the parameter in the D is transformed. ## calc : list of items to be calcuated ## 1: Calculate the conditional likelihood for subject i: f(Y_i|b_i). ## 2: Calculate the first derivative of item 1 above w.r.t $\beta$ and $\sigma$. ## 3: return the three parts of conditional mean's derivative w.r.t $\beta$ and $\sigma$. yi=datai$yi xi=datai$xi zi=datai\$zi obs<-!is.na(yi) xiobs<-as.matrix(xi[obs,,drop=F]) ziobs<-as.matrix(zi[obs,,drop=F]) yiobs<-as.vector(yi[obs]) niobs<-length(yiobs) if (transform==TRUE) D=exp(D) par = c(beta,D) etaiobs<-as.matrix(cbind(xiobs, t(t(ziobs)*qpt)))%*%as.matrix(par) ##print(t(t(ziobs)*qpt)) muiobs<-as.vector(exp(etaiobs)/(1+exp(etaiobs))) cond.lik<-numeric(niobs) ## conditional likelihood for each observation der.cond.lik<-numeric(length(par)) ## derivative of conditional likelihood. cond.lik = dbinom(yiobs, 1, muiobs) pbound = 10^(-10) cond.lik = ifelse(cond.lik100, 100, ifelse(y< -100, -100, y)) eta10=exp(eta10+gamma1*y) eta20=exp(eta20+gamma1*y) eta11=exp(eta11+gamma1*y) if( g==0 & gp==0) pg=1/(1+eta10+ eta20) else if (g==1 & gp==0) pg=eta10/(1+eta10+ eta20) else if (g==2 & gp==0) pg=eta20/(1+eta10+ eta20) else if (g==0 & gp==1) pg=1/(1+eta11) else if (g==1 & gp==1) pg=eta11/(1+eta11) else stop("There should be no such value of yp in the data") pg } ## code for some utility function fun.postrans <- function(u) { if (u) stop("Argument has to be positive") res <- log(u) res } fun.trans2pos <- function(u) { res<-exp(u) res } fun.rhotrans <- function(rho) { if (abs(rho) >1) stop("correlation coefficient must be between 1 and -1") tol= 1e-10 u = (1+rho)/2 if (u==0) u= u + tol if (u==1) u= 1- tol res<- log(u/(1-u)) res } fun.trans2rho <- function(u) { tol=10^308 if (u==Inf) u= tol if (u== -Inf) u= -tol rho<- 2*exp(u)/(1+exp(u))-1 rho }