inst/webtool/miscFunctions.R

In jcheng5/pqantimalarials: web tool for estimating under-five deaths caused by poor-quality antimalarials in sub-Saharan Africa

#library(rms)
#library(plyr)

getNormalLHS = function(N, mean, sd, aboveZeroOnly = TRUE) {
	# returns a LHS from the normal distibution with
	# mean and sd: splits the normal distribution
	# into N equiprobable sections and then returns a 
	# shuffled vector containing a sample from each section

	# N random numbers between 0 and 1/N
	randProb = runif(N)/N
	#print(randProb)
	
	# lower probability bounds for N equally weighted
	# probability sections
	sectionLowerBounds = 0:(N - 1)/N
	#print(sectionLowerBounds)
	
	# generates random probability taken from each
	# of the N equally weighted probability sections
	rand = randProb + sectionLowerBounds
	#print(rand)
	
	# reverse quantile function to get random samples
	# from each equally weighted section of the normal
	# probability distribution with mean and sd
	LHsample = qnorm(rand, mean, sd)
	#print("preRemoveNegatives")
	#print(LHsample)

	# if any values are below 0, set to 0 
	######## WARNING!!! : YOU MAY NOT WANT THIS!!! #############
	if (aboveZeroOnly) {
		lessThanZero <- LHsample > 0
		#print("lessThanZero")
		#print(lessThanZero)
	LHsample <- lessThanZero * LHsample
	}

	#print("postRemoveNegatives")
	#print(LHsample)

	# takes a sample from LHsample without replacement
	shuffledLHS = sample(LHsample)

	# return shuffled latin hypercube sample
	shuffledLHS
}

getUniformLHS = function(N, min, max) {
	# returns a LHS from the uniform distibution with
	# min and max: splits the uniform distribution
	# into N equiprobable sections and then returns a 
	# shuffled vector containing a sample from each section

	# N random numbers between 0 and 1/N
	randProb = runif(N)/N
	#print(randProb)
	
	# lower probability bounds for N equally weighted
	# probability sections
	sectionLowerBounds = 0:(N - 1)/N
	#print(sectionLowerBounds)
	
	# generates random probability taken from each
	# of the N equally weighted probability sections
	rand = randProb + sectionLowerBounds
	#print(rand)
	
	# reverse quantile function to get random samples
	# from each equally weighted section of the uniform
	# probability distribution with min and max
	LHsample = qunif(rand, min, max)
	#print(LHsample)
	
	# takes a sample from LHsample without replacement
	shuffledLHS = sample(LHsample)

	# return shuffled latin hypercube sample
	shuffledLHS
}

## Summary SE Summarizes data and is Taken from: http://www.cookbook-r.com/Manipulating_data/Summarizing_data/

## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data = NULL, measurevar, groupvars = NULL, na.rm = FALSE, conf.interval = 0.95, .drop = TRUE) {

	# New version of length which can handle NA's: if na.rm==T, don't count them
	length2 <- function(x, na.rm = FALSE) {
		if (na.rm) 
			sum(!is.na(x))
		else length(x)
	}

	# This is does the summary; it's not easy to understand...
	datac <- plyr::ddply(data, groupvars, .drop = .drop, 
					.fun = function(xx, col, na.rm) {
						c(N = length2(xx[, col], na.rm = na.rm),
						 mean = mean(xx[, col], na.rm = na.rm),
						 sd = sd(xx[, col], na.rm = na.rm),
						 max = max(xx[, col], na.rm = na.rm),
						 min = min(xx[, col], na.rm = na.rm),
						 q1 = quantile(xx[, col], na.rm = na.rm, 1/4),
						 q3 = quantile(xx[, col], na.rm = na.rm, 3/4),
						 median =  median(xx[, col], na.rm = na.rm))
					}, 
			measurevar, na.rm)

	# Rename the "mean" column    
	datac <- plyr::rename(datac, c("mean" = measurevar))

	datac$se <- datac$sd/sqrt(datac$N) # Calculate standard error of the mean

	# Confidence interval multiplier for standard error
	# Calculate t-statistic for confidence interval: 
	# e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
	ciMult <- qt(conf.interval/2 + 0.5, datac$N - 1)
	datac$ci <- datac$se * ciMult

	colnames(datac)[7] <- "q1"
	colnames(datac)[8] <- "q3"
	
	return(datac)
}

### Model ###
model = function(sales, fakePercent, deathRate) {
	#fakeProportion = fakePercent/100
	fakeProportion = fakePercent
	deaths = sales * fakeProportion * deathRate

	# add sum to end of vector
	deaths = c(deaths, sum(deaths))

	# return deaths
	deaths
}

# returns the standard deviation of a uniform distribution
# with a given min and max
stdUniform = function(min, max){
	std <- sqrt((1/12)*(max-min)^2)
}

counterfeitPRCC = function(x, sort.results = FALSE, sort.abs = FALSE) {
	# This code was adopted from Eili Klein <klein@cddep.org>
	
	# N = number of simulations
	N = length(x[, 1])
	#print("N:")
	#print(N)
	
	# k = number of input parameters + 1 output var
	k = length(x[1, ])
	#print("K:")
	#print(k)

	# Rank each input parameter
	r = matrix(NA, nrow = N, ncol = k)
	for (i in 1:k) {
		#r[, i] = rank(ties.method = "min",x[, i])
		r[, i] = rank(x[, i])
	}

	# save output ranks
	outputvar = r[, k]
	
	# r is the matrix where each column (1 for each input parameter) 
	# contains ranks of the simulation values for that parameter
	r = r[, -k]
	k = k - 1

	# If two of the input parameters have exactly the same ranking for
	# every run, then only one of the parameters should be used in the
	# calculation of PRCC
	dropcols = 0
	for (i in 1:(k - 1)) {
		dups = seq(1, k)
		for (j in 1:i) dups[j] = 0
		for (j in (i + 1):k) {
			a = which(r[, i] == r[, j])
			if (length(a) != N) {
				dups[j] = 0
			}
		}
		
		# keep track of duplicate columns
		for (j in 1:k) {
			if (dups[j] > 0) {
				dropcols = c(dropcols, j)
			}
		}
	}
	
	# remove 0 as a dropcol	
	dropcols = dropcols[-1]
	
	# get unique column numbers
	dropcols = unique(dropcols)
	
	# sort dropcol list
	dropcols = sort(dropcols)
	
	# reverse list of columns to be dropped so that you
	# can drop them using the # as an index and it wont
	# interfere with subsequent drops
	dropcols = rev(dropcols)
	
	# drop duplicate columns
	if (length(dropcols) > 0) {
		for (i in 1:length(dropcols)) {
			r = r[, -dropcols[i]]
		}
	}
	
	# adjust K if you droped any columns
	k = length(r[1, ])

	# bind output rankings to input rankings matrix
	r = cbind(r, outputvar)
	
	### Generate C Matrix ###
	mu = (1 + N)/2
	C.ij = matrix(NA, nrow = k + 1, ncol = k + 1)
	for (i in 1:(k + 1)) {
		for (j in 1:(k + 1)) {
			C.ij[i, j] = sum((r[, i] - mu) * (r[, j] - mu))/sqrt(sum((r[, i] - 
				mu)^2) * sum((r[, j] - mu)^2))
		}
	}

	B = rms::matinv(C.ij)

	gamma.ij = rep(0, k)
	t.iy = rep(0, k)
	p.iy = rep(0, k)
	
	for (i in 1:(k)) {
		# the PRCC between the ith input parameter and the outcome variable
		gamma.ij[i] = -B[i, k + 1]/sqrt(B[i, i] * B[k + 1, k + 1])
		
		# the significance of a nonzero PRCC is tested by computing t.iy
		# the distribution of t.iy approximates a students T with N-2 degrees of freedom
		t.iy[i] = gamma.ij[i] * sqrt((N - 2)/(1 - gamma.ij[i]))
		
		p.iy[i] = 2 * pt(-abs(t.iy[i]), df = N - 2)
	}

	# account for dropped columns
	dropcols = sort(dropcols)
	if (length(dropcols) > 0) {
		for (i in 1:length(dropcols)) {
			if (dropcols[i] > length(gamma.ij)) {
				# append to end
				gamma.ij = c(gamma.ij[1:(dropcols[i] - 1)], 0)
				t.iy = c(t.iy[1:(dropcols[i] - 1)], "dropped")
				p.iy = c(p.iy[1:(dropcols[i] - 1)], "-")
			} else {
				# insert in location
				gamma.ij = c(gamma.ij[1:(dropcols[i] - 1)], 0, gamma.ij[dropcols[i]:(length(gamma.ij))])
				t.iy = c(t.iy[1:(dropcols[i] - 1)], "dropped", t.iy[(dropcols[i]):(length(t.iy))])
				p.iy = c(p.iy[1:(dropcols[i] - 1)], "-", p.iy[(dropcols[i]):(length(p.iy))])
			}
		}
	}
	
	# calculate absolute value column to be used for sorting
	gamma.ij.abs = abs(gamma.ij)
	
	# create output dataframe
	vals = data.frame(cbind(gamma.ij, t.iy, p.iy, gamma.ij.abs), row.names = names(x[1:(length(x[1, 
		]) - 1)]))

	if (sort.results == TRUE) {
		if (sort.abs == TRUE) {
			vals = vals[order(vals$gamma.ij.abs, decreasing = TRUE), ]
		} else {
			vals = vals[order(vals$gamma.ij, decreasing = TRUE), ]
		}
	}

	# drop absolute value column
	vals = vals[, -4]
	
	vals

}