R/mm_model.R

In nicolemunoz99/multisensory: Predicting the conditions under which multisensory integration is favored

#' Optimal Integration Plot
#'
#' Plots "favorability of integration" (sensu Munoz & Blumstein, "Optimal Integration")
#' examining user-specified parameters. See Munoz & Blumstein's "Optimal Integration"
#' (under review at Behavioral Ecology), for derivation and list of assumptions.
#'
#' @param xlist A list having structure and named components equal to "optint.empty". See
#'   vignette for further requirements of xlist.opt .
#'   COL is Parameter varied x-axis of plot(s).
#'   ROW is Parameter varied between plots.
#'   SER is parameter varied within a plot.
#'   DIM1 is parameter of dimension-1 of parameter area for calculating favorability of integration.
#'   DIM2 is parameter of dimension-2 of parameter area for calculating favorability of integration.
#'   CONST1 is parameter held constant.
#'   CONST2 is parameter held constant.
#'   CONST3 is parameter held constant.
#'
#' @return Plot(s) of favorability of integration having user specified dimensions and values. The
#'   number of plots returned is equal to length(xlist$ROW$value)
#'
#' @examples
#' oiplot(optint.ex1)
#' oiplot(optint.ex2)
#'
#' @export


oiplot <- function(xlist){
# plots "favorability of integration" as a function of 8 parameters (3 of which, are held constant)
# as specified by user. Requires a list having the structure and named components like "optint.empty"

	if (!is.list(xlist)) stop('xlist must be a list with structure and named components like optint.empty')

	############################################################
	# associate parameter with plot feature from user input
	############################################################
	COL <- xlist$COL$param 		# name of parameter varied along x-axis
	ROW <- xlist$ROW$param 		# name of parameter varied across rows of graphs
	SER <- xlist$SER$param 		# name of parameter varied within a graph (i.e., different lines)
	DIM1 <- xlist$DIM1$param 	# name of parameter of one dimension of A matrix
	DIM2 <- xlist$DIM2$param 	# name of parameter of the other dimension of A matrix
	CONST1 <- xlist$CONST1$param	# name of parameter held constant
	CONST2 <- xlist$CONST2$param	# name of parameter held constant
	CONST3 <- xlist$CONST3$param	# name of parameter held constant

	if (length(COL) > 80) stop ('Maximum length for COL is 80')
	if (length(ROW) > 3) stop ('Maximum length for ROW is 3')
	if (length(SER) > 6) stop ('Maximum length for ROW is 6')
	if (length(DIM1) > 100 || length(DIM2) > 100) ("Maximum length for DIM1, DIM2 is 100")
	if (length(CONST1) != 1 || length(CONST2) != 1 || length(CONST3) != 1) stop('CONST1, CONST2,
	  CONST3 must have lengths = 1')

	# Are all parameters included in xlist?
	param.nam = c("bb.b", "u1", "u2", "prior.b", "bb.a", "k1", "k2", "s1")
	tt = unlist(lapply(xlist, function(x) x[[1]]))
	if (any(is.na(match(tt, param.nam)))) stop('xlist must contain name of parameters as in
	  optint.empty (i.e., u1, u2, prior.b, s1, k1, k2, bb.a, bb.b)')

	# Are all parameter values in xlist numeric?
	ss = lapply(xlist, function(x) x[[2]])
	nonnum = sapply(ss, function(x) any(is.numeric(x)))
	if (any(nonnum==F)) {
		exceptions = tt[which(nonnum==F)]
		stop(paste(exceptions, ' has been specified with non-numeric values'))
	}

	############################################################
	# define parameters based on list from user input
	############################################################
	param.name <- lapply(xlist, function(x){
		x$param
	})

	# prior probability of world in state b
	prior.b <- unlist(xlist[[which(param.name=="prior.b")]]["value"])


	# fraction of overlap of b and a distributions; assume 10% < overlap < %90
	u1 <- unlist(xlist[[which(param.name=="u1")]]["value"])
	u2 <- unlist(xlist[[which(param.name=="u2")]]["value"])


	# benefit of correct behavior in state "b" (bb.b) and correct behavior in state "a" (bb.a)
	bb.b <- unlist(xlist[[which(param.name=="bb.b")]]["value"])
	bb.a <- unlist(xlist[[which(param.name=="bb.a")]]["value"])

	# cost of using first and second stimulus
	k1 <- unlist(xlist[[which(param.name=="k1")]]["value"])
	k2 <- unlist(xlist[[which(param.name=="k2")]]["value"])

	# magnitude of first stimulus
	s1 <- unlist(xlist[[which(param.name=="s1")]]["value"])

	if(any(prior.b <= 0) || any (prior.b >= 1)) stop ('prior.b must be between 0 and 1')
	if(any(u1 <= 0) || any (u1 >= 1)) stop ('u1 must be between 0 and 1')
	if(any(u2 <= 0) || any (u2 >= 1)) stop ('u2 must be between 0 and 1')
	if(any(bb.b < 0)) stop ('bb.b must be >= 0')
	if(any(bb.a < 0)) stop ('bb.a must be >= 0')
	if(k1 < 0) stop ('k1 must be >= 0')
	if(k2 < 0) stop ('k2 must be >= 0')
	if(any(s1 < -10) || any(s1 > 10) || any(bb.b > 10) || any(bb.a > 10) || any(k1 > 10) || any(k2 > 10))
		warning ('Sensitivity of Favorability of Integration is best visualized for: -10 < s1 < 10;
		0 < bb.a < 10; 0 < bb.b < 10; 0 < k1 < 10; 0 < k2 < 10')


	############################################################
	# create "param" a table of parameter sets for calculating Value of Information
	############################################################


	col <- eval(parse(text = COL))
	row <- eval(parse(text = ROW))
	ser <- eval(parse(text = SER))
	xx <- eval(parse(text = DIM1))
	yy <- eval(parse(text = DIM2))
	const1 <- eval(parse(text = CONST1))
	const2 <- eval(parse(text = CONST2))
	const3 <- eval(parse(text = CONST3))

	# x-axis values
	x1 <- sapply(col, function(pp){
		rep(pp, length(row)*length(ser)*length(xx)*length(yy))
	})
	x1.a <- array(x1)

	# values of row-varying parameter
	x1.sub1 <- sapply(row, function(rrr){
		rep(rrr, length(ser)*length(xx)*length(yy))
	})
	x1.sub1.a <- rep(array(x1.sub1), length(col))

	# values of series-varying parameter
	x1.sub2 <- sapply(ser, function(ss){
		rep(ss, length(xx)*length(yy))
	})
	x1.sub2.a <- rep(unlist(x1.sub2), length(col)*length(row))


	# values of dimension 1 of matrix-varying parameter
	x1.sub3.a <- rep(xx, length(col)*length(row)*length(ser)*length(yy))

	# values of dimension 2 of matrix-varying parameter
	x1.sub4 <- sapply(yy, function(yyy){rep(yyy, length(xx))})
	x1.sub4.a <- rep(array(x1.sub4), length(row)*length(col)*length(ser))


	param1 <- cbind(x1.a, x1.sub1.a, x1.sub2.a, x1.sub3.a, x1.sub4.a)
	colnames(param1) <- cn <- sapply(c(COL, ROW, SER, DIM1, DIM2), function(x) x)

	# add values of constants to "param1"
	param <- cbind(param1, rep(const1, nrow(param1)), rep(const2, nrow(param1)), rep(const3, nrow(param1)))
	colnames(param)[6:8] <- c(CONST1, CONST2, CONST3)


	# rownames(param) <- sapply(1:nrow(param), function(x){
	# 	str_c(sapply(cn, function(y){ paste(y, "=", round(param[x,y],2), sep="")}), collapse=" ")
	# })


	############################################################
	# apply function "voi" to calculate integration results
	# results stored in arrays: vi1.u, vi2.u, int1.u, int2.u, st.use.u
	############################################################

	vi2.u <- vi1.u <- int2.u <- int1.u <- st.use.u <- array(NA, nrow(param)) # initiate unlisted arrays

	p.list <- vector("list", ncol(param)) # list to hold parameters that will be passed to function "voi"
	names(p.list) <- paste("t.", colnames(param), sep="")


	for(gg in 1:nrow(param)){

		p.list[1:ncol(param)] = param[gg,] # pass parameter values to p.list

		temp = do.call("voi", p.list) # evaluate voi with parameters in p.list

		int1.u[gg] <- temp$int1[1] # is the 1st stimulus used?
		vi1.u[gg] <- temp$vi.1[1] # value of information of 1st stimulus
		int2.u[gg] <- temp$int2[1] # is the 2nd stimulus used?
		vi2.u[gg] <- temp$vi.2[1] # value of information of 2nd stimulus
		st.use.u[gg] <- temp$st.use[1] # 1--> only first stimulus used; 2 --> only 2nd stimulus used; 3 --> integration

		temp = NA # reset
	}



	############################################################
	# rearrange integration output into lists for plotting
	############################################################



	skel <- lapply(col, function(a){ # structure for relisting integration arrays back into matrices suitable for graphing
		x = lapply(row, function(b){
			y = lapply(ser, function(c){
				z = matrix(NA, length(xx), length(yy))
				rownames(z) = paste(DIM1, round(xx,2), sep = "=")
				colnames(z) = paste(DIM2, round(yy,2), sep = "=")
				z
			}); names(y) = paste(SER, ser, sep = "="); y
		}); names(x) = paste(ROW, round(row,2), sep = "="); x
	})

	names(skel) <- paste(COL, round(col,2), sep = "=")

	int2 <- relist(int2.u, skel) # execute relisting
	int1 <- relist(int1.u, skel)
	vi2 <- relist(vi2.u, skel)
	vi1 <- relist(vi1.u, skel)
	st.use <- relist(st.use.u, skel)




	############################################################
	# make plots
	############################################################


	ul.st.use <- unlist(unlist(st.use, recursive = F), recursive = F)

	for(ii in c(3)){ # only look at integration (i.e., values in st.use that equal 3)

		# proportion of parameter space where a given stimulus/stimuli is/are used
		prop = sapply(ul.st.use, function(x){
			length(which(x == ii))/length(x)
		})
		# relist prop with structure similar to st.use
		skel2 = lapply(col, function(x){ # skeleton for relisting prop similar to int2
			t2 = lapply(row, function(y){
				t1 = lapply(ser, function(z){
					NA
				}); names(t1) = paste(SER, round(ser,2), sep="="); t1
			}); names(t2) = paste(ROW, round(row,2), sep="="); t2
		}); names(skel2) = paste(COL, round(col,2), sep="=")
		areas = relist(prop, skel2) # relist prop.3

		# restructure areas list in structure suitable for plotting
		y.areas <- lapply(1:length(ser), function(x){
			t2 = lapply(1:length(row), function(y){
				t1 = sapply(1:length(col), function(z){
					areas[[z]][[y]][[x]]
				}); names(t1) = paste(COL, round(col,2), sep="="); t1
			}); names(t2) = paste(ROW, round(row,2), sep="="); t2
		}); names(y.areas) = paste(SER, round(ser,2), sep="=")


		par(mfrow = c(1,1), cex.axis = 1, mar = c(6, 5, 1, 8), las = 1, lwd = 1, xpd = T)

		# graph stuff
		for(RR in 1:length(row)){
			# plot each COL in a new plot
			plot(col, y.areas[[1]][[RR]], xlab = COL, ylab = "A (favorability of integration)", ylim = c(0,1),  type = "l", cex.lab = 1)
			if(length(ser) > 1){
				lapply(2:length(ser), function(x){
					points(col, y.areas[[x]][[RR]], type = "l", lty=x)
				})
			}
			# add legend
			leg1 = legend(1.06*max(col), 0.9, legend = round(ser,3), title = SER, lty = 1:length(ser), cex = 1)
			# add labels
			text(leg1$rect$left, 0.3, paste(ROW, row[RR], sep=" = "), adj = 0, cex = 1.2)
			text(leg1$rect$left, 0.0, paste(CONST1, const1, sep=" = "), adj = 0, cex = 1.0)
			text(leg1$rect$left, 0.1, paste(CONST2, const2, sep=" = "), adj = 0, cex = 1.0)
			text(leg1$rect$left, 0.2, paste(CONST3, const3, sep=" = "), adj = 0, cex = 1.0)

			# text(leg1$rect$left, 1, paste("st.use =", ii), adj = 0, cex = 1.5)
		} # end plot "for" loop for a given outcome (RR)
	} # end loop over all possible outcomes (ii)

}