R/PH_samplegamma.r

###########################################################################################
# This code is used to estimate the parameters of a given data set using the MRH
#	methodology.  One to many covariates can be included under the proportional 
#	hazards assumption.
###########################################################################################

PH_samplegamma = function(Mval, delta, X, Xstdz, outfilename, burnIn, maxIter, thin, k,
GR, fix.burnIn, fix.thin, fix.max, writenum, sysseconds, systemtime, checknum, n, betaNames,
numParams, betaLB, betaUB, beta, betas.LB.stdz, betas.UB.stdz, stdIndex, Xmeans, Xsds, mat01,
RmpInd, one_RmpInd, inBin, failBin, mvec, Rmp, a, lambda, RmpNoPruneIndx, RmpRowNums, initialValues,
continue.chain, gamma.init, loopctr.start){

    # Initialize and set parameters
	c = d = 1
	
	#########################################################
	# Start output file for parameter estimates
	#########################################################
	# Write out ds, betas, H00, Rmps
    if(continue.chain == FALSE){
        gamma.mp = rep(0.5, 2^Mval-1)
        if(GR == TRUE){ gamma.mp[RmpNoPruneIndx] = runif(length(RmpNoPruneIndx), .2, .8)    }
       if(numParams > 0){
            if(length(betaNames) != numParams){	betaNames = 1:numParams	}
                write.table(matrix(c('iteration', paste('d',1:2^Mval, sep = ''), paste('beta', betaNames, sep = '.'),
                    'H00', paste('Rmp', RmpRowNums, sep = ''), paste('gammamp', RmpRowNums, sep = ''), 'a', 'lambda'), nrow = 1),
                    paste(outfilename, '/MCMCchains.txt', sep = ''), row.names = FALSE, col.names = FALSE)
        } else {
            write.table(matrix(c('iteration', paste('d',1:2^Mval, sep = ''),
					'H00', paste('Rmp', RmpRowNums, sep = ''), paste('gammamp', RmpRowNums, sep = ''), 'a', 'lambda'), nrow = 1),
					paste(outfilename, '/MCMCchains.txt', sep = ''), row.names = FALSE, col.names = FALSE)
        }
    } else {
        gamma.mp = gamma.init
    }
    initialValues = c(initialValues, list(gamma.init = gamma.mp))
    # number of rounding values needed: 2^M (d) + length(betaNames) (optional) + 1 (H00) + 2^M-1 (Rmp) +
    # 1 (a) + 1 (lambda) + 2^M-1 (gamma)
    round.values = rep(16, 2^Mval + 1 + 2^Mval-1 + 1 + 1 + 2^Mval-1)
    if(numParams > 0){  round.values = c(round.values, rep(16, length(betaNames)))   }

	#########################################################
	# Run Gibbs sampler
	#########################################################
	convergence = best.thin.found = FALSE
	if(fix.thin == TRUE){ best.thin.found = TRUE }
	outdata = NULL
	loopctr = loopctr.start
	
	# assessIndex keeps the index of the parameters that need to be tested for convergence and plotted
	# is equal to the betas, 1 H, 2^M-1 Rmps minus those that are not pruned, 2^M-1 gammas 
	assessIndex = c(1:(numParams+1), (1:(2^Mval-1))[RmpNoPruneIndx]+numParams+1, (1:(2^Mval-1))[RmpNoPruneIndx] + 2^Mval-1+numParams+1)+2^Mval+1
	
	while((loopctr <= maxIter & convergence == FALSE & fix.max == FALSE) | 
		  (loopctr <= maxIter & fix.max == TRUE)){
    
		######## Draw H00 ########
		# Calculate F, which is needed in the posterior of H00
		F = calc_F(mat01 = mat01, Rmp = Rmp, M = Mval, inBinMat = inBin)
		# Draw H00
		H00 = rgamma(1, shape = a+sum(delta), scale = (1/lambda+sum(exp(X%*%beta)*F))^-1)
	
		######## Draw Rmp values ########
		# Calculate H, which is needed in the posterior of the Rmps
		H = calc_H(mat01 = mat01, H00 = H00, Rmp = Rmp, M = Mval, inBinMat = inBin)
		# Draw Rmps
		for(Rmpctr in RmpNoPruneIndx){
			Rmp[Rmpctr] = .Call("arms", c(.001, .999), f = function(x) logRmpPost_PHcovs(x, RmpFull = Rmp, H00val = H00, kval = k, 
																						 aval = a, 
																						 gamma.mpval = gamma.mp[Rmpctr],
																						 betavec = beta, X.mat = X, 
																						 mval = mvec[Rmpctr], 
																						 RmpIndic = RmpInd[,Rmpctr], 
																						 one_RmpIndic = one_RmpInd[,Rmpctr],
																						 deltavals = delta, Mvalue = Mval, 
																						 inBinMat = inBin, mat01 = mat01, 
																						 formula = Rmpctr), Rmp[Rmpctr], 
								as.integer(1), new.env())
		}
		
		##### Draw beta ########
		# H needs to be recalculated first
		if(numParams > 0){
			if(length(stdIndex) > 0){
				betas.stdz = beta*Xsds
				H00.stdz = H00*exp(sum(betas.stdz/Xsds*Xmeans))
				H = calc_H(mat01 = mat01, H00 = H00.stdz, Rmp = Rmp, M = Mval, inBinMat = inBin)
				for(betaCtr in 1:numParams){
					H00.stdz = H00*exp(sum(betas.stdz/Xsds*Xmeans))
					H = calc_H(mat01 = mat01, H00 = H00.stdz, Rmp = Rmp, M = Mval, inBinMat = inBin)
					betas.stdz[betaCtr] = .Call("arms", c(betas.LB.stdz[betaCtr], betas.UB.stdz[betaCtr]), 
													   f = function(x) logbetaPost_PH(x, betaFull = betas.stdz, 
																					  deltavals = delta, Xmatrix = Xstdz, 
																					  Hvals = H, whichBeta = betaCtr, mu.beta = 0, 
																					  sigma.beta = 10), betas.stdz[betaCtr], 
													   as.integer(1), new.env())
				}
				beta = betas.stdz/Xsds
			} else {
				H = calc_H(mat01 = mat01, H00 = H00, Rmp = Rmp, M = Mval, inBinMat = inBin)
				for(betaCtr in 1:numParams){
					beta[betaCtr] = .Call("arms", c(betaLB[betaCtr], betaUB[betaCtr]), 
										  f = function(x) logbetaPost_PH(x, betaFull = beta, deltavals = delta, 
																		 Xmatrix = X, Hvals = H, whichBeta = betaCtr, mu.beta = 0, 
																		 sigma.beta = 10), beta[betaCtr], 
										  as.integer(1), new.env())
				}
			}
		}
		
		##### Draw gamma #########
		for(gammactr in RmpNoPruneIndx){
			gamma.mp[gammactr] = .Call("arms", c(.001, .999), 
									   f = function(x) loggammaPost(x, Rmp = Rmp[gammactr], kval = k, aval = a,
																	cval = c, dval = d, mvecval = mvec[gammactr]),
									   gamma.mp[gammactr], as.integer(1), new.env())
		}
		
		##### Draw lambda ########
		lambda = .Call("arms", c(.0001, 10), f = function(x) logLambdaPost(x, mu.lval = 100, H00val = H00, aval = a),
					   lambda, as.integer(1), new.env())
		##### Draw a #########
		a = sample.a(1:50, lambdaval = lambda, mu.aval = 4, H00val = H00, kval = k, 
							 mvec = mvec, Mval = Mval, Rmpvec = Rmp, gamma.mpvec = gamma.mp)

		#######################################################################################
		#					Gather and provide information for user 
		#######################################################################################
		# At 100 iterations, give the user a run time estimate
        if(loopctr == loopctr.start + 99){ checkRunTime(maxIter-loopctr.start+1, systemtime)  }
		if((loopctr %% thin) == 0 & loopctr >= burnIn){
			d.iter = calc_dfast(Mval, Rmp, H00, mat01)
			if(numParams > 0){
				outdata = rbind(outdata, matrix(c(loopctr, d.iter, beta, H00, Rmp, gamma.mp, a, lambda), nrow = 1))
			} else {
				outdata = rbind(outdata, matrix(c(loopctr, d.iter, H00, Rmp, gamma.mp, a, lambda), nrow = 1))
			}
		}
        # If the routine reaches 100,000 iterations, check that the rounding values for the MCMC
        # text file will work (not too large)
        if(loopctr == 100000){
            testData = read.table(paste(outfilename, '/MCMCchains.txt', sep = ''), header = TRUE)
            round.values = FindRoundVals(testData)
        }
        # Every 10,000 or writenum iterations, write out the data set and empty the working one.
        if((loopctr %% writenum) == 0 & loopctr >= burnIn){
            if(nrow(outdata) > 1){
                outdata = cbind(outdata[,1], sapply(2:ncol(outdata), function(x) round(outdata[,x], round.values[x-1])))
            } else {
                outdata = matrix(c(outdata[,1], sapply(2:ncol(outdata), function(x) round(outdata[,x], round.values[x-1]))), nrow = 1)
            }
            write.table(outdata, paste(outfilename, '/MCMCchains.txt', sep = ''),
            col.names = FALSE, row.names = FALSE, append = TRUE)
            if(loopctr != maxIter){	outdata = NULL	}
        }
        # Every 100,000 or checknum iterations, check for convergence with the full data set (minus the burned values)
        if((loopctr %% checknum) == 0 & loopctr >= burnIn){
            chainRes = checkMCMCParams(best.thin.found, assessIndex, thin, burnIn, round.values, outfilename,
            convergence, fix.burnIn)
            thin = chainRes$thin
            burnIn = chainRes$burnIn
            convergence = chainRes$convergence
            best.thin.found = chainRes$best.thin.found
        }
		if((loopctr %% 5000) == 0){	print(paste(loopctr, 'iterations completed...'))	}
		
		# Increment loopctr
		loopctr = loopctr+1
	}
	outputdata = list(assessIndex = assessIndex, numberofEstParams = 2 + 1 + length(RmpNoPruneIndx)*2 + numParams,
                        numParams = numParams, failBin = failBin, inBin = inBin, GR = GR, loopctr = loopctr,
                        convergence = convergence, initialValues = initialValues, burnIn = burnIn, thin = thin)
	return(outputdata)
	
}

Try the MRH package in your browser

Any scripts or data that you put into this service are public.

MRH documentation built on May 2, 2019, 11:10 a.m.