mutost.calibrate: Calibrating 'mutost' for ABC
In olli0601/abc.star: abc.star

Description Usage Arguments Value Note References See Also Examples

Calibrate the one-sample equivalence test for population means of normal summary values for ABC inference. The one-sample mutost can be used within ABC to test the null hypothesis that the underlying population mean of the simulated summary values is not similar to the sample mean of the observed summary values. It is applicable when the simulated and observed summary values follow a normal distribution, or when normality cannot be rejected.

Different types of calibrations are available, see Notes for details:

(what=ALPHA) compute the ABC false positive rate for given critical region,
(what=CR) calibrate the critical region for given ABC false positive rate,
(what=MXPW) calibrate the critical region and the equivalence region for given ABC false positive rate and maximum power,
(what=KL) calibrate the critical region, the equivalence region and the number of simulated summary values for given ABC false positive rate, maximum power and sample standard deviation of the observed data.

In the ideal case, the calibration KL is used. However, the KL calibration requires multiple i. i. d. instances of observed summary statistics at each ABC iteration. If this is not available, the MXPW calibration should be used.

Depending on the type of calibration, some of the following inputs must be specified (see Examples).

mutost.calibrate(n.of.x = NA, s.of.x = NA, n.of.y = n.of.x, s.of.y = NA,
  what = "MXPW", c.u = NA, tau.u = NA, tau.u.ub = NA, mx.pw = 0.9,
  alpha = 0.01, max.it = 100, pow_scale = 1.5, tol = 1e-05,
  debug = FALSE, plot = FALSE, plot_debug = FALSE, verbose = FALSE)

`n.of.x`	Number of observed summary values
`s.of.x`	Standard deviation of observed summary values
`n.of.y`	Number of simulated summary values
`s.of.y`	Standard deviation of simulated summary values
`what`	Character string to indicate the type of calibration to be performed
`c.u`	Upper boundary point of the critical region
`tau.u`	Upper boundary point of the equivalence region
`tau.u.ub`	Guess on the upper boundary point of the equivalence region
`mx.pw`	Maximum power at the point of equality
`alpha`	Level of the equivalence test
`max.it`	Maximum number of optimization steps at each calibration hierarchy
`pow_scale`	Scale for the support of the standardized power. The power is truncated to `pow_scale*[-tau.u,tau.u]` and then standardized
`tol`	Required error tolerance in calibrating the actual maximum power to the requested maximum power
`debug`	Flag to switch off C implementation
`plot`	Flag to plot calibrations
`plot_debug`	Flag to plot at each calibration iteration
`verbose`	Flag to run in verbose mode

vector

(what=ALPHA) This calibration requires the inputs c.u, tau.u with c.u<tau.u and c.u>0. The output contains the corresponding ABC false positive rate alpha. This option does not specify any of the free ABC parameters, but may be useful to determine the ABC false positive rate for uncalibrated ABC routines.
(what=CR) This calibration requires the input tau.u, alpha with tau.u>0 and default alpha=0.01. The output contains the corresponding critical region [c.l, c.u], which corresponds to the ABC tolerance region typically denoted by [-epsilon, epsilon]. The resulting critical region may be empty under this test, if stochasticity is too large. This is an intermediate calibration step and may result in unsuitable power properties (see Examples).
(what=MXPW) This calibration requires the inputs alpha, mx.pw, with default values 0.01 and 0.9 respectively. The output contains the corresponding critical region [c.l, c.u] (to be used in ABC, see Notes on (2)), and the corresponding equivalence region [tau.l, tau.u] that gives a suitable ABC accept/reject probability if the simulated summary values are close to the observed summary values. As a check to the numerical calibrations, the actual power at the point of equality is returned (pw.cmx).
(what=KL) This calibration can be used when a set of observed summary values is available. It is desirable because it specifies the number of simulated summary values so that the power is very close to the desired summary likelihood in terms of the KL divergence. The inputs are alpha, mx.pw, with default values 0.01 and 0.9 respectively. values so that the power is very close to the desired summary likelihood in terms of the KL divergence. Thus, the output consists of the corresponding critical region [c.l, c.u] (to be used in ABC, see Notes on (2)), the equivalence region [tau.l, tau.u], and the number of simulated summary values needed (n.of.y). As a check to the numerical calibrations, the KL divergence is returned (KL). It is desirable to compare the power to the summary likelihood in terms of the KL divergence, see References.

Note that the underlying test statistic depends on the sample standard deviation of the simulated summary values, s.of.y. Consequently, the free ABC parameters must be re-calibrated when s.of.y changes.

The lower boundary point of the critical region c.l is always fixed to -c.u, because this choice implies that the power is maximized at the point of equality rho=0. This is also commonly used in uncalibrated ABC routines.

http://arxiv.org/abs/1305.4283

vartest.calibrate, ztest.calibrate, ratetest.calibrate

yn			<- 60
ymean 		<- 1.3
ysigma 		<- 0.8
sim 		<- rnorm(yn, ymean, ysigma)
sim 		<- (sim - mean(sim))/sd(sim) * ysigma + ymean

# Example 1A: calibrate critical region and power of ABC accept/reject step (default)
# this requires to specify alpha, mx.pw (calibration parameters), 
# n.of.y, s.of.y (summary parameters), 
# tau.u.ub (for numerical optimization)
# note: this is the default calibration because it specifies all ABC parameters
# for sensible calibration parameters, and only requires minimal information on 
# the observed summary values. If a set of observed summary values is available,
# use the calibrations in Example 3.

tau.u.ub 	<- 0.5
ans 		<- mutost.calibrate(n.of.y=length(sim), s.of.y=sd(sim), 
	tau.u.ub=tau.u.ub, what='MXPW', mx.pw=0.9, alpha=0.01, plot=TRUE)
						
# Example 1B: calibrate critical region and power of ABC accept/reject step
# note: critical region depends on ysigma

ysigma 		<- 1.3
sim 		<- rnorm(yn, ymean, ysigma)
sim 		<- (sim - mean(sim))/sd(sim) * ysigma + ymean
ans 		<- mutost.calibrate(n.of.y=length(sim), s.of.y=sd(sim), what='MXPW', mx.pw=0.9, alpha=0.01, plot=TRUE)
						
# Example 2A: calculate ABC false positive rate for given ABC tolerance
# this requires to specify c.u, tau.u (ad-hoc ABC parameters), n.of.y, s.of.y (summary parameters)
# note: can be useful to compute the ABC false positive rate for uncalibrated ABC routines

cus			<- c(0.1, 0.2, 0.3, 0.4, 0.5)
tau.u		<- 0.6817	
ans 		<- sapply(cus, function(c.u)	mutost.calibrate(n.of.y=length(sim), s.of.y=sd(sim), c.u=c.u, tau.u=tau.u, what='ALPHA')	) 

# Example 2B: calibrate critical region for given ABC false positive rate and equivalence region
# this requires to specify alpha (calibration parameters), tau.u (ad-hoc ABC parameter), n.of.y, s.of.y (summary parameters)
# note: this is just an intermediate calibration and may result in unsuitable power properties

ans 		<- mutost.calibrate(n.of.y=length(sim), s.of.y=sd(sim), tau.u=1.2,
	what='CR', alpha=0.01, plot=TRUE)

# Example 3A: calibrate critical region, power of ABC accept/reject step, and #simulated data points
# this requires to specify alpha, mx.pw (calibration parameters), 
# n.of.x, s.of.x, s.of.y (summary parameters), 
# tau.u.ub (for numerical optimization)
# note: the advantage here is that the KL divergence is also minimised, but we also
# need to have a set of observed summary values
#
# case sd(sim)>sd(obs): the power function is not "too tight" for yn=xn.

xn 			<- 60
xmean 		<- 1
xsigma 		<- 1
obs 		<- rnorm(xn, xmean, xsigma)
obs 		<- (obs - mean(obs))/sd(obs) * xsigma + xmean
ans 		<- mutost.calibrate(n.of.x=length(obs), s.of.x= sd(obs), s.of.y=sd(sim), 
	what='KL', mx.pw=0.9, alpha=0.01, plot=TRUE, debug=TRUE)

# case sd(sim)<sd(obs): the power function is "too tight" for yn=xn.

ysigma 		<- 0.5
sim 		<- rnorm(yn, ymean, ysigma)
sim 		<- (sim - mean(sim))/sd(sim) * ysigma + ymean
ans 		<- mutost.calibrate(n.of.x=length(obs), s.of.x= sd(obs), s.of.y=sd(sim), 
	what='KL', mx.pw=0.9, alpha=0.01, plot=TRUE, debug=TRUE)

system.time({ for(i in 1:1e3) mutost.calibrate(n.of.x=length(obs), s.of.x= sd(obs), s.of.y=sd(sim), what='KL', mx.pw=0.9, alpha=0.01) 	})

## Not run: 
abc.presim.uprior.mu<- function(abc.nit, xn, xmean, xsigma, prior.l, prior.u, ysigma, yn=NA )		
{		
	ans			<- vector("list",5)
	names(ans)	<- c("x","xn","xmean","xsigma","sim")
	obs 		<- rnorm(xn, xmean, xsigma)
	obs 		<- (obs - mean(obs))/sd(obs) * xsigma + xmean
	ans[["x"]]			<- obs
	ans[["xmean"]]		<- xmean
	ans[["xsigma"]]		<- xsigma
	
	ans[["sim"]]		<- sapply(1:abc.nit, function(i)
			{					
				ymu		<- runif(1, prior.l, prior.u)
				y		<- rnorm(yn, ymu, sd=ysigma)
				tmp		<- c(yn, ymu, ysigma, mean(y), sd(y) )									
				tmp					
			})								
	rownames(ans[["sim"]])	<- c('m','ymu','ysigma','ysmean','yssd')
	ans
}


abc.presim	<- abc.presim.uprior.mu( 	abc.nit=1e7, xn=60, xmean=1.34, xsigma=1.4, 
		prior.l=1.34-5, prior.u=1.34+5, ysigma=1.4, yn=60 )
abc.df		<- as.data.table(t(abc.presim$sim))								
abc.df[, it:=seq_len(nrow(abc.df))]								
tmp			<- abc.df[,  as.list( mutost.calibrate(	n.of.y=m, s.of.y=yssd, tau.u.ub=1, what='MXPW', mx.pw=0.9, alpha=0.01)[1:4] ), by='it']
abc.df		<- merge(abc.df, tmp, by='it')								
abc.df[, T:= abc.df[, ysmean-1.34]]
abc.df[, ABC_C90:= abc.df[, c.l<=T & T<=c.u]]

ggplot( subset(abc.df, ABC_C90), aes(x=ymu-abc.presim$xmean)) + geom_histogram()
ggplot( subset(abc.df, ABC_C90), aes(x=ymu-abc.presim$xmean)) + geom_histogram(aes(y= ..density..))

## End(Not run)