R/abc.R
In bcm-uga/abc: Tools for Approximate Bayesian Computation (ABC)
######################################################################
#
# abc.R
#
# copyright (c) 2011-05-30, Katalin Csillery, Olivier Francois and
# Michael GB Blum with some initial code from Mark Beaumont
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License,
# version 3, as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but without any warranty; without even the implied warranty of
# merchantability or fitness for a particular purpose. See the GNU
# General Public License, version 3, for more details.
#
# A copy of the GNU General Public License, version 3, is available
# at http://www.r-project.org/Licenses/GPL-3
#
# Part of the R/abc package
# Contains: abc, return.abc, print.abc, is.abc, summary.abc, hist.abc, plot.abc
#
######################################################################
abc <- function(target, param, sumstat, tol, method, hcorr = TRUE,
transf = "none", logit.bounds = c(0,0), subset = NULL,
kernel = "epanechnikov", numnet = 10, sizenet = 5,
lambda = c(0.0001,0.001,0.01), trace = FALSE, maxit =
500, ...){
call <- match.call()
## general checks that the function is used correctly
## ###################################################
if(missing(target)) stop("'target' is missing")
if(missing(param)) stop("'param' is missing")
if(missing(sumstat)) stop("'sumstat' is missing")
if(!is.matrix(param) && !is.data.frame(param) && !is.vector(param)) stop("'param' has to be a matrix, data.frame or vector.")
if(!is.matrix(sumstat) && !is.data.frame(sumstat) && !is.vector(sumstat)) stop("'sumstat' has to be a matrix, data.frame or vector.")
if(missing(tol)) stop("'tol' is missing")
if(missing(method)) stop("'method' is missing with no default")
if(!any(method == c("rejection", "loclinear", "neuralnet","ridge"))){
stop("Method must be rejection, loclinear, or neuralnet or ridge")
}
if(method == "rejection") rejmethod <- TRUE
else rejmethod <- FALSE
if(!any(kernel == c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine"))){
kernel <- "epanechnikov"
warning("Kernel is incorrectly defined. Setting to default (Epanechnikov)")
}
if(is.data.frame(param)) param <- as.matrix(param)
if(is.data.frame(sumstat)) sumstat <- as.matrix(sumstat)
if(is.list(target)) target <- unlist(target)
if(is.vector(sumstat)) sumstat <- matrix(sumstat, ncol=1)
if(length(target)!=dim(sumstat)[2]) stop("Number of summary statistics in 'target' has to be the same as in 'sumstat'.")
## stop if zero var in sumstat
## #########################
nss <- length(sumstat[1,])
cond1 <- !any(as.logical(apply(sumstat, 2, function(x) length(unique(x))-1)))
if(cond1) stop("Zero variance in the summary statistics.")
## transformations
## ################
ltransf <- length(transf)
if(is.vector(param)){
numparam <- 1
param <- matrix(param, ncol=1)
}
else numparam <- dim(param)[2]
for (i in 1:ltransf){
if(sum(transf[i] == c("none","log","logit")) == 0){
stop("Transformations must be none, log, or logit.")
}
if(transf[i]=="logit"){
if(logit.bounds[i,1] >= logit.bounds[i,2]){
stop("Logit bounds are incorrect.")
}
}
}
## if no logit change logit.bounds to NULL
## if(any(transf) == "logit") logit.bounds <- NULL
## no transformation should be applied when rejmethod is true
if(rejmethod){
if(!all(transf == "none")){
warning("No transformation is applied when the simple rejection is used.", call.=F)
}
transf[1:numparam] <- "none"
}
else{
if(numparam != ltransf){
if(length(transf) == 1){
transf <- rep(transf[1], numparam)
warning("All parameters are \"", transf[1], "\" transformed.", sep="", call.=F)
}
else stop("Number of parameters is not the same as number of transformations.", sep="", call.=F)
}
}
## parameter and/or sumstat values that are to be excluded
## #######################################################
gwt <- rep(TRUE,length(sumstat[,1]))
gwt[attributes(na.omit(sumstat))$na.action] <- FALSE
if(missing(subset)) subset <- rep(TRUE,length(sumstat[,1]))
gwt <- as.logical(gwt*subset)
## extract names of parameters and statistics if given
## ###################################################
if(!length(colnames(param))){
warning("No parameter names are given, using P1, P2, ...")
paramnames <- paste("P", 1:numparam, sep="")
}
else paramnames <- colnames(param)
if(!length(colnames(sumstat))){
warning("No summary statistics names are given, using S1, S2, ...")
statnames <- paste("S", 1:nss, sep="")
}
else statnames <- colnames(sumstat)
## scale everything
## #################
scaled.sumstat <- sumstat
for(j in 1:nss){
scaled.sumstat[,j] <- normalise(sumstat[,j],sumstat[,j][gwt])
}
for(j in 1:nss){
target[j] <- normalise(target[j],sumstat[,j][gwt])
}
## calculate euclidean distance
## ############################
sum1 <- 0
for(j in 1:nss){
sum1 <- sum1 + (scaled.sumstat[,j]-target[j])^2
}
dist <- sqrt(sum1)
## includes the effect of gwt in the tolerance
dist[!gwt] <- floor(max(dist[gwt])+10)
## wt1 defines the region we're interested in
ceiling(length(dist)*tol)->nacc
sort(dist)[nacc]->ds
wt1 <- (dist <= ds)
aux<-cumsum(wt1)
wt1 <- wt1 & (aux<=nacc)
if(kernel == "gaussian")
{
wt1 <- rep(TRUE, length(dist))
}
## transform parameters
## ######################
for (i in 1:numparam){
if(transf[i] == "log"){
if(min(param[,i]) <= 0){
cat("log transform: values out of bounds - correcting...")
x.tmp <- ifelse(param[,i] <= 0,max(param[,i]),param[,i])
x.tmp.min <- min(x.tmp)
param[,i] <- ifelse(param[,i] <= 0, x.tmp.min,param[,i])
}
param[,i] <- log(param[,i])
}
else if(transf[i] == "logit"){
if(min(param[,i]) <= logit.bounds[i,1]){
x.tmp <- ifelse(param[,i] <= logit.bounds[i,1],max(param[,i]),param[,i])
x.tmp.min <- min(x.tmp)
param[,i] <- ifelse(param[,i] <= logit.bounds[i,1], x.tmp.min,param[,i])
}
if(max(param[,i]) >= logit.bounds[i,2]){
x.tmp <- ifelse(param[,i] >= logit.bounds[i,2],min(param[,i]),param[,i])
x.tmp.max <- max(x.tmp)
param[,i] <- ifelse(param[,i] >= logit.bounds[i,2], x.tmp.max,param[,i])
}
param[,i] <- (param[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
param[,i] <- log(param[,i]/(1-param[,i]))
}
} # end of parameter transformations
## select summary statistics in region
## ###################################
ss <- sumstat[wt1,]
unadj.values <- param[wt1,]
## if simple rejection or in the selected region there is no var in sumstat
## ########################################################################
statvar <- as.logical(apply(cbind(sumstat[wt1,]), 2, function(x) length(unique(x))-1))
cond2 <- !any(statvar)
if(cond2 && !rejmethod)
stop("Zero variance in the summary statistics in the selected region. Try: checking summary statistics, choosing larger tolerance, or rejection method.")
if(rejmethod){
if(cond2) warning("Zero variance in the summary statistics in the selected region. Check summary statistics, consider larger tolerance.")
weights <- rep(1,length=sum(wt1))
adj.values <- NULL
residuals <- NULL
lambda <- NULL
}
else{
if(cond2) cat("Warning messages:\nStatistic(s)", statnames[!statvar], "has/have zero variance in the selected region.\nConsider using larger tolerance or the rejection method or discard this/these statistics, which might solve the collinearity problem in 'lsfit'.\n", sep=", ")
## weights
if(kernel == "epanechnikov") weights <- 1 - (dist[wt1]/ds)^2
if(kernel == "rectangular") weights <- dist[wt1]/ds
if(kernel == "gaussian") weights <- 1/sqrt(2*pi)*exp(-0.5*(dist/(ds/2))^2)
if(kernel == "triangular") weights <- 1 - abs(dist[wt1]/ds)
if(kernel == "biweight") weights <- (1 - (dist[wt1]/ds)^2)^2
if(kernel == "cosine") weights <- cos(pi/2*dist[wt1]/ds)
## regression correction
## ######################
if(method == "loclinear"){
fit1 <- lsfit(scaled.sumstat[wt1,],param[wt1,],wt=weights)
#pred <- t(fit1$coeff) %*% c(1,target)
pred <- t(structure(cbind(fit1$coefficients)[fit1$qr$pivot,], names=names(fit1$coefficients))) %*% c(1,target)
pred <- matrix(pred, ncol=numparam, nrow=sum(wt1), byrow=TRUE)
#residuals <- fit1$residuals
residuals<-param[wt1,]-t(t(structure(cbind(fit1$coefficients)[fit1$qr$pivot,], names=names(fit1$coefficients))) %*%t(cbind(1,scaled.sumstat[wt1,])))
residuals<-cbind(residuals)
#rrr<<-residuals
the_m<-apply(residuals,FUN=mean,2)
residuals<-sapply(1:numparam,FUN=function(x){residuals[,x]-the_m[x]})
pred <- sapply(1:numparam,FUN=function(x){pred[,x]+the_m[x]})
sigma2<-apply(as.matrix(residuals),FUN=function(x){sum((x)^2 *weights)/sum(weights)},MARGIN=2)
aic<-sum(wt1)*sum(log(sigma2))+2*(nss+1)*numparam
bic<-sum(wt1)*sum(log(sigma2))+log(sum(wt1))*(nss+1)*numparam
if(hcorr == TRUE){
fit2 <- lsfit(scaled.sumstat[wt1,],log(residuals^2),wt=weights)
auxaux<-t(structure(cbind(fit2$coefficients)[fit2$qr$pivot,], names=names(fit2$coefficients))) %*% c(1,target)
pred.sd <- sqrt(exp(auxaux))
pred.sd <- matrix(pred.sd, nrow=sum(wt1), ncol=numparam, byrow=T)
pred.si<-t(t(structure(cbind(fit2$coefficients)[fit2$qr$pivot,], names=names(fit2$coefficients))) %*%t(cbind(1,scaled.sumstat[wt1,])))
pred.si<-sqrt(exp(pred.si))
#residuals2<-log(residuals^2)-t(t(structure(cbind(fit2$coefficients)[fit2$qr$pivot,], names=names(fit2$coefficients))) %*%t(cbind(1,scaled.sumstat[wt1,])))
#fv <- sqrt(exp(log(fit1$residuals^2) - fit2$residuals))
#fv <- sqrt(exp(residuals2))
# adj.values <- pred + (pred.sd*residuals) / fv # correction heteroscedasticity
adj.values <- pred + (pred.sd*residuals) / pred.si
residuals<-(pred.sd*residuals) / pred.si
}
else{
adj.values <- pred + residuals
}
colnames(adj.values) <- colnames(unadj.values)
lambda <- NULL
}
## neural network regression
## ##########################
if(method == "neuralnet"){
linout <- TRUE
## normalise parameters
param.mad <- c()
for(i in 1:numparam){
param.mad[i] <- mad(param[,i][gwt]) # save for renormalisation
param[,i] <- normalise(param[,i],param[,i][gwt])
}
lambda <- sample(lambda, numnet, replace=T)
fv <- array(dim=c(sum(wt1), numparam, numnet))
pred <- matrix(nrow=numparam, ncol=numnet)
for (i in 1:numnet){
fit1 <- nnet(scaled.sumstat[wt1,], param[wt1,], weights = weights, decay = lambda[i],
size = sizenet, trace = trace, linout = linout, maxit = maxit, ...)
cat(i)
fv[,,i] <- fit1$fitted.values
pred[,i] <- predict(fit1, data.frame(rbind(target)))
}
cat("\n")
pred.med <- apply(pred, 1, median)
pred.med <- matrix(pred.med, nrow=sum(wt1), ncol=numparam, byrow=T)
fitted.values <- apply(fv, c(1,2), median)
residuals <- param[wt1,] - fitted.values # median of fitted values par nnets for each accepted point and parameter
if(hcorr == TRUE){
pred2 <- matrix(nrow=numparam, ncol=numnet)
fv2 <- array(dim=c(sum(wt1), numparam, numnet))
for (i in 1:numnet){
fit2 <- nnet(scaled.sumstat[wt1,], log(residuals^2), weights = weights, decay = lambda[i], size = sizenet, trace = trace, linout = linout, ...)
cat(i)
fv2[,,i] <- fit2$fitted.values
pred2[,i] <- predict(fit2, data.frame(rbind(target)))
}
cat("\n")
pred.sd <- sqrt(exp(apply(pred2, 1, median)))
pred.sd <- matrix(pred.sd, nrow=sum(wt1), ncol=numparam, byrow=T)
fv.sd <- sqrt(exp(apply(fv2, c(1,2), median)))
adj.values <- pred.med + (pred.sd*residuals)/fv.sd # correction heteroscedasticity
residuals<-(pred.sd*residuals)/fv.sd
}
else{
adj.values <- pred.med + residuals
}
colnames(adj.values) <- colnames(unadj.values)
## renormalise
for(i in 1:numparam){
## residuals[,i] <- residuals[,i]*param.mad[i] not much sense...
adj.values[,i] <- adj.values[,i]*param.mad[i]
}
} # end of neuralnet
if(method == "ridge"){
#print("la")
## normalise parameters
param.mad <- c()
for(i in 1:numparam){
param.mad[i] <- mad(param[,i][gwt]) # save for renormalisation
param[,i] <- normalise(param[,i],param[,i][gwt])
}
# print("ici")
numnet<-length(lambda)
#lambda <- sample(lambda, numnet, replace=T)
fv <- array(dim=c(sum(wt1), numparam, numnet))
pred <- matrix(nrow=numparam, ncol=numnet)
mataux<-sqrt(diag(weights))
paramaux<-as.matrix(mataux%*%param[wt1,])
#print(dim(paramaux))
scaledaux<-mataux%*%scaled.sumstat[wt1,]
#targetaux<-drop(mataux%*%target
for(parcur in (1:numparam))
{
#cat(parcur)
fit1<-lm.ridge(paramaux[,parcur]~scaledaux,lambda=lambda)
for (i in 1:numnet){
# fit1 <- nnet(scaled.sumstat[wt1,], param[wt1,], weights = weights, decay = lambda[i],
# size = sizenet, trace = trace, linout = linout, maxit = maxit, ...)
# cat(i)
fv[,parcur,i] <- drop(cbind(1,scaled.sumstat[wt1,])%*%(rbind(coef(fit1))[i,]))
pred[parcur,i] <- drop(c(1, target) %*%(rbind(coef(fit1))[i,]))
}
}
# cat("\n")
pred.med <- apply(pred, 1, median)
pred.med <- matrix(pred.med, nrow=sum(wt1), ncol=numparam, byrow=T)
fitted.values <- apply(fv, c(1,2), median)
residuals <- param[wt1,] - fitted.values # median of fitted values par nnets for each accepted point and parameter
if(hcorr == TRUE){
pred2 <- matrix(nrow=numparam, ncol=numnet)
fv2 <- array(dim=c(sum(wt1), numparam, numnet))
for(parcur in (1:numparam))
{
lresidaux<-(mataux%*%(log(residuals[,parcur]^2)))
fit2<-lm.ridge(lresidaux~scaledaux,lambda=lambda)
for (i in 1:numnet){
#cat(i)
fv2[,parcur,i] <- drop(cbind(1,scaled.sumstat[wt1,])%*%(rbind(coef(fit2))[i,]))
pred2[parcur,i] <- drop(c(1, target) %*%(rbind(coef(fit2))[i,]))
}
}
cat("\n")
pred.sd <- sqrt(exp(apply(pred2, 1, median)))
pred.sd <- matrix(pred.sd, nrow=sum(wt1), ncol=numparam, byrow=T)
fv.sd <- sqrt(exp(apply(fv2, c(1,2), median)))
adj.values <- pred.med + (pred.sd*residuals)/fv.sd # correction heteroscedasticity
residuals<- (pred.sd*residuals)/fv.sd
}
else{
adj.values <- pred.med + residuals
}
colnames(adj.values) <- colnames(unadj.values)
## renormalise
for(i in 1:numparam){
## residuals[,i] <- residuals[,i]*param.mad[i] not much sense...
adj.values[,i] <- adj.values[,i]*param.mad[i]
}
} # end of ridge
} # end of else than rejmethod
## back transform parameter values
## ################################
if(numparam == 1){
unadj.values <- matrix(unadj.values, ncol=1)
if(method!="rejection")
{adj.values <- matrix(adj.values, ncol=1)
residuals <- matrix(residuals, ncol=1)}
}
for (i in 1:numparam){
if(transf[i] == "log"){
unadj.values[,i] <- exp(unadj.values[,i])
adj.values[,i] <- exp(adj.values[,i])
}
else if(transf[i] == "logit"){
unadj.values[,i] <- exp(unadj.values[,i])/(1+exp(unadj.values[,i]))
unadj.values[,i] <- unadj.values[,i]*(logit.bounds[i,2]-logit.bounds[i,1])+logit.bounds[i,1]
adj.values[,i] <- exp(adj.values[,i])/(1+exp(adj.values[,i]))
adj.values[,i] <- adj.values[,i]*(logit.bounds[i,2]-logit.bounds[i,1])+logit.bounds[i,1]
}
}
abc.return(transf, logit.bounds, method, call, numparam, nss, paramnames, statnames,
unadj.values, adj.values, ss, weights, residuals, dist, wt1, gwt, lambda, hcorr,aic,bic)
}
abc.return <- function(transf, logit.bounds, method, call, numparam, nss, paramnames, statnames,
unadj.values, adj.values, ss, weights, residuals, dist, wt1, gwt, lambda, hcorr,aic,bic){
if(method == "rejection"){
out <- list(unadj.values=unadj.values, ss=ss, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="rejection", numparam=numparam, numstat=nss, names=list(parameter.names=paramnames, statistics.names=statnames))
}
else if(method == "loclinear"){
out <- list(adj.values=adj.values, unadj.values=unadj.values,
ss=ss, weights=weights, residuals=residuals, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="loclinear", hcorr = hcorr, numparam=numparam, numstat=nss,aic=aic,bic=bic,
names=list(parameter.names=paramnames, statistics.names=statnames))
}
else if(method == "ridge"){
out <- list(adj.values=adj.values, unadj.values=unadj.values,
ss=ss, weights=weights, residuals=residuals, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="ridge", hcorr = hcorr, numparam=numparam, numstat=nss,
names=list(parameter.names=paramnames, statistics.names=statnames))
}
else if(method == "neuralnet"){
out <- list(adj.values=adj.values, unadj.values=unadj.values,
ss=ss, weights=weights, residuals=residuals, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="neuralnet", hcorr = hcorr, lambda=lambda, numparam=numparam, numstat=nss,
names=list(parameter.names=paramnames, statistics.names=statnames))
}
class(out) <- "abc"
out
}
is.abc <- function(x){
if (inherits(x, "abc")) TRUE
else FALSE
}
print.abc <- function(x, ...){
if (!inherits(x, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- x
cl <- abc.out$call
cat("Call:\n")
dput(cl, control=NULL)
cat("Method:\n")
## rejection
if (abc.out$method == "rejection")
cat("Rejection\n\n")
## loclinear
else if (abc.out$method == "loclinear"){
if(abc.out$hcorr == T){
cat("Local linear regression\n")
cat("with correction for heteroscedasticity\n\n")
}
else cat("Local linear regression\n\n")
}
## ridge
else if (abc.out$method == "ridge"){
if(abc.out$hcorr == T){
cat("Ridge regression\n")
cat("with correction for heteroscedasticity\n\n")
}
else cat("Ridge regression\n\n")
}
## nnet
else if (abc.out$method == "neuralnet"){
if(abc.out$hcorr == T){
cat("Non-linear regression via neural networks\n")
cat("with correction for heteroscedasticity\n\n")
}
else cat("Non-linear regression via neural networks\n\n")
}
cat("Parameters:\n")
cat(abc.out$names$parameter.names, sep=", ")
cat("\n\n")
cat("Statistics:\n")
cat(abc.out$names$statistics.names, sep=", ")
cat("\n\n")
cat("Total number of simulations", length(abc.out$na.action), "\n\n")
cat("Number of accepted simulations: ", dim(abc.out$unadj.values)[1], "\n\n")
invisible(abc.out)
}
summary.abc <- function(object, unadj = FALSE, intvl = .95, print = TRUE, digits = max(3, getOption("digits")-3), ...){
if (!inherits(object, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- object
np <- abc.out$numparam
cl <- abc.out$call
parnames <- abc.out$names[[1]]
npri <- dim(abc.out$unadj.values)[1]
npost <- dim(abc.out$unadj.values)[1]
myprobs <- c((1-intvl)/2, 0.5, 1-(1-intvl)/2)
if(print){
cat("Call: \n")
dput(cl, control=NULL)
}
if(abc.out$method == "rejection" || unadj){
if(print) cat(paste("Data:\n abc.out$unadj.values (",npost," posterior samples)\n\n", sep=""))
res <- matrix(abc.out$unadj.values, ncol=np)
mymin <- apply(res, 2, min)
mymax <- apply(res, 2, max)
mymean <- apply(res, 2, mean)
mymode <- apply(res, 2, getmode, ...)
quants <- apply(res, 2, quantile, probs = myprobs)
sums <- rbind(mymin, quants[1,], quants[2,], mymean, mymode, quants[3,], mymax)
dimnames(sums) <- list(c("Min.:", paste(c(myprobs[1])*100, "% Perc.:", sep=""),
"Median:", "Mean:", "Mode:",
paste(c(myprobs[3])*100, "% Perc.:", sep=""), "Max.:"),
parnames)
}
else{
if(print){
cat(paste("Data:\n abc.out$adj.values (",npost," posterior samples)\n", sep=""))
cat(paste("Weights:\n abc.out$weights\n\n", sep=""))
}
res <- matrix(abc.out$adj.values, ncol=np)
wt <- abc.out$weights
mymin <- apply(res, 2, min)
mymax <- apply(res, 2, max)
mymean <- apply(res, 2, weighted.mean, w = wt)
mymode <- apply(res, 2, getmode, wt/sum(wt), ...)
quants <- apply(res, 2, function(x) rq(x~1, tau = myprobs, weights = wt)$coef)
sums <- rbind(mymin, quants[1,], quants[2,], mymean, mymode, quants[3,], mymax)
dimnames(sums) <- list(c("Min.:", paste("Weighted ", c(myprobs[1])*100, " % Perc.:", sep=""),
"Weighted Median:", "Weighted Mean:", "Weighted Mode:",
paste("Weighted ", c(myprobs[3])*100, " % Perc.:", sep=""), "Max.:"),
parnames)
}
class(sums) <- "table"
if(print){
if(length(digits) == 1)
print(round(sums, digits = digits), quote=FALSE)
else
print(apply(rbind(digits, sums), 2, function(a) round(a[-1], digits = a[1])), quote=FALSE)
invisible(sums)
}
else sums
}
getmode <- function(x, weights = NULL, ...){
## if(missing(bw) | missing(kernel) | missing(window) | missing(n)) warning("density.default() was used to calculate the posterior mode")
d <- density(x, weights = weights)
d$x[which(d$y == max(d$y))][1]
}
hist.abc <- function(x, unadj = FALSE, true = NULL,
file = NULL, postscript = FALSE, onefile = TRUE, ask = !is.null(deviceIsInteractive()),
col.hist = "grey", col.true = "red", caption = NULL, ...){
if (!inherits(x, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- x
np <- abc.out$numparam
parnames <- abc.out$names[[1]]
if(is.null(caption)) mycaption <- as.graphicsAnnot(parnames)
else mycaption <- caption
## checks if true is given
if(!is.null(true)){
if(!is.vector(true)){
stop("Supply true parameter value(s) as a vector.", call.=F)
}
if (length(true) != np){
stop("Number of true values has to be the same as number of parameters in 'x'.", call.=F)
}
cond <- isTRUE(c(match(names(true), parnames), match(names(true), parnames)))
if(cond) stop("Names do not match in 'x' and 'true'.", call.=F)
}
if(abc.out$method == "rejection") res <- abc.out$unadj.values
else if (unadj){
res <- abc.out$unadj.values
}
else res <- abc.out$adj.values
## Devices
## ##########
save.devAskNewPage <- devAskNewPage()
if(!is.null(file)){
file <- substitute(file)
if(!postscript) pdf(file = paste(file, "pdf", sep="."), onefile=onefile)
if(postscript) postscript(file = paste(file, "ps", sep="."), onefile=onefile)
}
else{
if (ask && prod(par("mfcol")) < np) {
devAskNewPage(TRUE)
}
}
myhist <- list()
for(i in 1:np){
myxlim <- range(res[,i], true[i])
myhist[[i]] <- hist(res[,i], xlab = parnames[i], main = paste("Posterior histogram of ", mycaption[i], sep=""),
col = col.hist, xlim = myxlim, ...)
myhist[[i]]$xname <- parnames[i]
if(!is.null(true)){
abline(v = true[i], lwd = 2, col = col.true)
}
}
if(!is.null(file)){
dev.off()
}
else devAskNewPage(save.devAskNewPage)
names(myhist) <- parnames
invisible(myhist)
} # end of hist.abc
##########################################################################################
## plots: for each parameter
## 1. prior density
## 2. posterior + prior density
## 3. distances vs parameters
## 4. histogram of the residuals
##########################################################################################
plot.abc <- function(x, param, subsample = 1000, true = NULL,
file = NULL, postscript = FALSE, onefile = TRUE, ask = !is.null(deviceIsInteractive()), ...){
if (!inherits(x, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- x
mymethod <- abc.out$method
if(mymethod == "rejection")
stop("Diagnostic plots can be displayed only when method is \"loclinear\", \"neuralnet\" or \"ridge\".", cal.=F)
if(!is.matrix(param) && !is.data.frame(param) && !is.vector(param)) stop("'param' has to be a matrix, data.frame or vector.", call.=F)
if(is.null(dim(param))) param <- matrix(param, ncol=1)
if(is.data.frame(param)) param <- as.matrix(param)
np <- abc.out$numparam
numsim <- length(param)/np
alldist <- log(abc.out$dist)
myregion <- abc.out$region
residuals <- abc.out$residuals
parnames <- abc.out$names$parameter.names
transf <- abc.out$transf
logit.bounds <- abc.out$logit.bounds
##check if param is compatible with x
cond <- isTRUE(c(match(colnames(param), parnames), match(parnames, colnames(param))))
if(cond) stop("'abc.out' and 'param' are not compatible; paramater names are different.", call.=F)
## checks if true is given
if(!is.null(true)){
if(!is.vector(true)){
stop("Supply true parameter value(s) as a vector.", call.=F)
}
if (length(true) != np){
stop("Number of true values has to be the same as number of parameters in 'x'.", call.=F)
}
cond <- isTRUE(c(match(names(true), parnames), match(names(true), parnames)))
if(cond) stop("Names do not match in 'x' and 'true'.", call.=F)
}
rej <- abc.out$unadj.values
res <- abc.out$adj.values
if(np == 1) rej <- matrix(rej, ncol=1)
if(is.vector(param)){
np.orig <- 1
nsim <- length(param)
}
else if(is.matrix(param)){
np.orig <- dim(param)[2]
nsim <- dim(param)[1]
myorder <- match(parnames, colnames(param))
if(isTRUE(myorder-1:np)){
param <- param[, myorder]
warning("'param' is being re-ordered according to 'abc.out'...", call.=F, immediate.=T)
}
}
## check if param has the right dimensions
if(np.orig != np){
stop("The number parameters supplied in \"param\" is different from that in \"x\".", call.=F)
}
## for (i in 1:np){
## if(transf[i] == "log"){
## if(min(param[,i]) <= 0){
## warning("Correcting out of bounds values for \"log\" transformed parameters.", call.=F, immediate.=T)
## x.tmp <- ifelse(param[,i] <= 0,max(param[,i]),param[,i])
## x.tmp.min <- min(x.tmp)
## param[,i] <- ifelse(param[,i] <= 0, x.tmp.min,param[,i])
## }
## param[,i] <- log(param[,i])
## if(min(rej[,i]) <= 0){
## warning("Correcting out of bounds values for \"log\" transformed parameters.", call.=F, immediate.=T)
## x.tmp <- ifelse(rej[,i] <= 0,max(rej[,i]),rej[,i])
## x.tmp.min <- min(x.tmp)
## rej[,i] <- ifelse(rej[,i] <= 0, x.tmp.min,rej[,i])
## }
## rej[,i] <- log(rej[,i])
## if(min(res[,i]) <= 0){
## warning("Correcting out of bounds values for \"log\" transformed parameters.", call.=F, immediate.=T)
## x.tmp <- ifelse(res[,i] <= 0,max(res[,i]),res[,i])
## x.tmp.min <- min(x.tmp)
## res[,i] <- ifelse(res[,i] <= 0, x.tmp.min,res[,i])
## }
## res[,i] <- log(res[,i])
## if(!is.null(true)){
## if(min(true[i]) <= 0){
## x.tmp <- ifelse(true[i] <= 0,max(true[i]),true[i])
## x.tmp.min <- min(x.tmp)
## true[,i] <- ifelse(true[i] <= 0, x.tmp.min,true[i])
## }
## true[i] <- log(true[i])
## }
## }
## else if(transf[i] == "logit"){
## if(min(param[,i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(param[,i] <= logit.bounds[i,1],max(param[,i]),param[,i])
## x.tmp.min <- min(x.tmp)
## param[,i] <- ifelse(param[,i] <= logit.bounds[i,1], x.tmp.min,param[,i])
## }
## if(max(param[,i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(param[,i] >= logit.bounds[i,2],min(param[,i]),param[,i])
## x.tmp.max <- max(x.tmp)
## param[,i] <- ifelse(param[,i] >= logit.bounds[i,2], x.tmp.max,param[,i])
## }
## param[,i] <- (param[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## param[,i] <- log(param[,i]/(1-param[,i]))
## if(min(rej[,i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(rej[,i] <= logit.bounds[i,1],max(rej[,i]),rej[,i])
## x.tmp.min <- min(x.tmp)
## rej[,i] <- ifelse(rej[,i] <= logit.bounds[i,1], x.tmp.min,rej[,i])
## }
## if(max(rej[,i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(rej[,i] >= logit.bounds[i,2],min(rej[,i]),rej[,i])
## x.tmp.max <- max(x.tmp)
## rej[,i] <- ifelse(rej[,i] >= logit.bounds[i,2], x.tmp.max,rej[,i])
## }
## rej[,i] <- (rej[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## rej[,i] <- log(rej[,i]/(1-rej[,i]))
## if(min(res[,i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(res[,i] <= logit.bounds[i,1],max(res[,i]),res[,i])
## x.tmp.min <- min(x.tmp)
## res[,i] <- ifelse(res[,i] <= logit.bounds[i,1], x.tmp.min,res[,i])
## }
## if(max(res[,i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(res[,i] >= logit.bounds[i,2],min(res[,i]),res[,i])
## x.tmp.max <- max(x.tmp)
## res[,i] <- ifelse(res[,i] >= logit.bounds[i,2], x.tmp.max,res[,i])
## }
## res[,i] <- (res[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## res[,i] <- log(res[,i]/(1-res[,i]))
## if(!is.null(true)){
## if(min(true[i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(true[i] <= logit.bounds[i,1],max(true[i]),true[i])
## x.tmp.min <- min(x.tmp)
## true[i] <- ifelse(true[i] <= logit.bounds[i,1], x.tmp.min,true[i])
## }
## if(max(true[i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(true[i] >= logit.bounds[i,2],min(true[i]),true[i])
## x.tmp.max <- max(x.tmp)
## true[i] <- ifelse(true[i] >= logit.bounds[i,2], x.tmp.max,true[i])
## }
## true[i] <- (true[i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## true[i] <- log(true[i]/(1-true[i]))
## }
## }
## } # end of parameter transformations
## devices
## ########
save.devAskNewPage <- devAskNewPage()
if(!is.null(file)){
file <- substitute(file)
if(!postscript) pdf(file = paste(file, "pdf", sep="."), onefile=onefile)
if(postscript) postscript(file = paste(file, "ps", sep="."), onefile=onefile)
}
else{
if (ask && 1 < np) {
devAskNewPage(TRUE)
}
}
## if param is too large, draw a random sample from it the size of
## which can be set by the user
mysample <- sample(1:nsim, subsample)
## plots
## #####
# prior, reg, rej, true
ltys <- c(3, 1, 1, 3)
lwds <- c(1, 2, 1, 2)
cols <- c("black", "red", "black", "black")
unadj <- TRUE
for(i in 1:np){
par(mfcol=c(2,2))
prior.d <- density(param[mysample,i])
post.d <- density(res[,i])
rej.d <- density(rej[,i])
myxlim <- range(c(post.d$x, rej.d$x, true))
myylim <- range(c(post.d$y, rej.d$y, prior.d$y))
par(cex = 1, cex.main = 1.2, cex.lab = 1.1, lwd = 2)
##if(transf[i] == "none")
myxlab <- parnames[i]
##else if(transf[i] == "log") myxlab <- paste("log(", parnames[i], ")", sep="")
##else if(transf[i] == "logit") myxlab <- paste("log(", parnames[i], "/(1-",parnames[i],"))", sep="")
## 1. prior
## ########
plot(prior.d, main = "", col=cols[1], lty = ltys[1], lwd=lwds[1], xlab=myxlab, ...)
title("Prior", sub = paste("N =", prior.d$n, " Bandwidth =", formatC(prior.d$bw)))
## 2. posterior
## ############
plot(post.d, col=cols[2], lty = ltys[2], lwd=lwds[2], xlim = myxlim, ylim = myylim, main = "", xlab=myxlab, ...)
lines(prior.d, col=cols[1], lty = ltys[1], lwd=lwds[1], ...)
lines(rej.d, col = cols[3], lty = ltys[3], lwd=lwds[3], ...)
if(!is.null(true)){
segments(x0 = true[i], x1 = true[i], y0 = myylim[1], y1 = myylim[2], col=cols[4], lty = ltys[4], lwd=lwds[4])
mtext(text="True value", side = 1, line=2, at=true[i])
}
title(paste("Posterior with \"", mymethod, "\"", sep=""), sub=paste("N =", post.d$n, " Bandwidth =", formatC(post.d$bw)))
mtext("\"rejection\" and prior as reference", side=3, line = 0.5)
## 3. distances
## ############
mypch <- 19
myylim <- range(alldist[mysample])#*c(1,1.2)
plot(param[mysample,i], alldist[mysample], col=cols[1],
xlab=myxlab, ylab="log Euclidean distance", ylim=myylim, pch=mypch, ...)
title("Euclidean distances")
mtext(paste("N(All / plotted) = ", numsim, "/", subsample, sep=" "), side = 3, line = .5)
points(param[myregion,i], alldist[myregion], col=cols[2], pch=mypch, ...)
mypar <- par()
if(!is.null(true)){
segments(x0 = true[i], x1 = true[i], y0 = myylim[1], y1 = myylim[2], col=cols[4], lty = ltys[4], lwd=lwds[4])
mtext(text="True value", side = 1, line=2, at=true[i])
}
## 4. residuals
## ############
if(mymethod=="loclinear") mymain <- "Residuals from lsfit()"
if(mymethod=="neuralnet") mymain <- "Residuals from nnet()"
if(mymethod=="ridge") mymain <- "Residuals from lm.ridge()"
qqnorm(residuals, pch=mypch, main=mymain, sub="Normal Q-Q plot", xlab="Theoretical quantiles", ylab="Residuals",...)
qqline(residuals)
} # np
par(mfcol=c(1,1))
if(!is.null(file)){
dev.off()
}
else devAskNewPage(save.devAskNewPage)
invisible()
} # end of plot.abc
bcm-uga/abc documentation built on May 11, 2019, 3:03 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
######################################################################
#
# abc.R
#
# copyright (c) 2011-05-30, Katalin Csillery, Olivier Francois and
# Michael GB Blum with some initial code from Mark Beaumont
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License,
# version 3, as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but without any warranty; without even the implied warranty of
# merchantability or fitness for a particular purpose. See the GNU
# General Public License, version 3, for more details.
#
# A copy of the GNU General Public License, version 3, is available
# at http://www.r-project.org/Licenses/GPL-3
#
# Part of the R/abc package
# Contains: abc, return.abc, print.abc, is.abc, summary.abc, hist.abc, plot.abc
#
######################################################################
abc <- function(target, param, sumstat, tol, method, hcorr = TRUE,
transf = "none", logit.bounds = c(0,0), subset = NULL,
kernel = "epanechnikov", numnet = 10, sizenet = 5,
lambda = c(0.0001,0.001,0.01), trace = FALSE, maxit =
500, ...){
call <- match.call()
## general checks that the function is used correctly
## ###################################################
if(missing(target)) stop("'target' is missing")
if(missing(param)) stop("'param' is missing")
if(missing(sumstat)) stop("'sumstat' is missing")
if(!is.matrix(param) && !is.data.frame(param) && !is.vector(param)) stop("'param' has to be a matrix, data.frame or vector.")
if(!is.matrix(sumstat) && !is.data.frame(sumstat) && !is.vector(sumstat)) stop("'sumstat' has to be a matrix, data.frame or vector.")
if(missing(tol)) stop("'tol' is missing")
if(missing(method)) stop("'method' is missing with no default")
if(!any(method == c("rejection", "loclinear", "neuralnet","ridge"))){
stop("Method must be rejection, loclinear, or neuralnet or ridge")
}
if(method == "rejection") rejmethod <- TRUE
else rejmethod <- FALSE
if(!any(kernel == c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine"))){
kernel <- "epanechnikov"
warning("Kernel is incorrectly defined. Setting to default (Epanechnikov)")
}
if(is.data.frame(param)) param <- as.matrix(param)
if(is.data.frame(sumstat)) sumstat <- as.matrix(sumstat)
if(is.list(target)) target <- unlist(target)
if(is.vector(sumstat)) sumstat <- matrix(sumstat, ncol=1)
if(length(target)!=dim(sumstat)[2]) stop("Number of summary statistics in 'target' has to be the same as in 'sumstat'.")
## stop if zero var in sumstat
## #########################
nss <- length(sumstat[1,])
cond1 <- !any(as.logical(apply(sumstat, 2, function(x) length(unique(x))-1)))
if(cond1) stop("Zero variance in the summary statistics.")
## transformations
## ################
ltransf <- length(transf)
if(is.vector(param)){
numparam <- 1
param <- matrix(param, ncol=1)
}
else numparam <- dim(param)[2]
for (i in 1:ltransf){
if(sum(transf[i] == c("none","log","logit")) == 0){
stop("Transformations must be none, log, or logit.")
}
if(transf[i]=="logit"){
if(logit.bounds[i,1] >= logit.bounds[i,2]){
stop("Logit bounds are incorrect.")
}
}
}
## if no logit change logit.bounds to NULL
## if(any(transf) == "logit") logit.bounds <- NULL
## no transformation should be applied when rejmethod is true
if(rejmethod){
if(!all(transf == "none")){
warning("No transformation is applied when the simple rejection is used.", call.=F)
}
transf[1:numparam] <- "none"
}
else{
if(numparam != ltransf){
if(length(transf) == 1){
transf <- rep(transf[1], numparam)
warning("All parameters are \"", transf[1], "\" transformed.", sep="", call.=F)
}
else stop("Number of parameters is not the same as number of transformations.", sep="", call.=F)
}
}
## parameter and/or sumstat values that are to be excluded
## #######################################################
gwt <- rep(TRUE,length(sumstat[,1]))
gwt[attributes(na.omit(sumstat))$na.action] <- FALSE
if(missing(subset)) subset <- rep(TRUE,length(sumstat[,1]))
gwt <- as.logical(gwt*subset)
## extract names of parameters and statistics if given
## ###################################################
if(!length(colnames(param))){
warning("No parameter names are given, using P1, P2, ...")
paramnames <- paste("P", 1:numparam, sep="")
}
else paramnames <- colnames(param)
if(!length(colnames(sumstat))){
warning("No summary statistics names are given, using S1, S2, ...")
statnames <- paste("S", 1:nss, sep="")
}
else statnames <- colnames(sumstat)
## scale everything
## #################
scaled.sumstat <- sumstat
for(j in 1:nss){
scaled.sumstat[,j] <- normalise(sumstat[,j],sumstat[,j][gwt])
}
for(j in 1:nss){
target[j] <- normalise(target[j],sumstat[,j][gwt])
}
## calculate euclidean distance
## ############################
sum1 <- 0
for(j in 1:nss){
sum1 <- sum1 + (scaled.sumstat[,j]-target[j])^2
}
dist <- sqrt(sum1)
## includes the effect of gwt in the tolerance
dist[!gwt] <- floor(max(dist[gwt])+10)
## wt1 defines the region we're interested in
ceiling(length(dist)*tol)->nacc
sort(dist)[nacc]->ds
wt1 <- (dist <= ds)
aux<-cumsum(wt1)
wt1 <- wt1 & (aux<=nacc)
if(kernel == "gaussian")
{
wt1 <- rep(TRUE, length(dist))
}
## transform parameters
## ######################
for (i in 1:numparam){
if(transf[i] == "log"){
if(min(param[,i]) <= 0){
cat("log transform: values out of bounds - correcting...")
x.tmp <- ifelse(param[,i] <= 0,max(param[,i]),param[,i])
x.tmp.min <- min(x.tmp)
param[,i] <- ifelse(param[,i] <= 0, x.tmp.min,param[,i])
}
param[,i] <- log(param[,i])
}
else if(transf[i] == "logit"){
if(min(param[,i]) <= logit.bounds[i,1]){
x.tmp <- ifelse(param[,i] <= logit.bounds[i,1],max(param[,i]),param[,i])
x.tmp.min <- min(x.tmp)
param[,i] <- ifelse(param[,i] <= logit.bounds[i,1], x.tmp.min,param[,i])
}
if(max(param[,i]) >= logit.bounds[i,2]){
x.tmp <- ifelse(param[,i] >= logit.bounds[i,2],min(param[,i]),param[,i])
x.tmp.max <- max(x.tmp)
param[,i] <- ifelse(param[,i] >= logit.bounds[i,2], x.tmp.max,param[,i])
}
param[,i] <- (param[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
param[,i] <- log(param[,i]/(1-param[,i]))
}
} # end of parameter transformations
## select summary statistics in region
## ###################################
ss <- sumstat[wt1,]
unadj.values <- param[wt1,]
## if simple rejection or in the selected region there is no var in sumstat
## ########################################################################
statvar <- as.logical(apply(cbind(sumstat[wt1,]), 2, function(x) length(unique(x))-1))
cond2 <- !any(statvar)
if(cond2 && !rejmethod)
stop("Zero variance in the summary statistics in the selected region. Try: checking summary statistics, choosing larger tolerance, or rejection method.")
if(rejmethod){
if(cond2) warning("Zero variance in the summary statistics in the selected region. Check summary statistics, consider larger tolerance.")
weights <- rep(1,length=sum(wt1))
adj.values <- NULL
residuals <- NULL
lambda <- NULL
}
else{
if(cond2) cat("Warning messages:\nStatistic(s)", statnames[!statvar], "has/have zero variance in the selected region.\nConsider using larger tolerance or the rejection method or discard this/these statistics, which might solve the collinearity problem in 'lsfit'.\n", sep=", ")
## weights
if(kernel == "epanechnikov") weights <- 1 - (dist[wt1]/ds)^2
if(kernel == "rectangular") weights <- dist[wt1]/ds
if(kernel == "gaussian") weights <- 1/sqrt(2*pi)*exp(-0.5*(dist/(ds/2))^2)
if(kernel == "triangular") weights <- 1 - abs(dist[wt1]/ds)
if(kernel == "biweight") weights <- (1 - (dist[wt1]/ds)^2)^2
if(kernel == "cosine") weights <- cos(pi/2*dist[wt1]/ds)
## regression correction
## ######################
if(method == "loclinear"){
fit1 <- lsfit(scaled.sumstat[wt1,],param[wt1,],wt=weights)
#pred <- t(fit1$coeff) %*% c(1,target)
pred <- t(structure(cbind(fit1$coefficients)[fit1$qr$pivot,], names=names(fit1$coefficients))) %*% c(1,target)
pred <- matrix(pred, ncol=numparam, nrow=sum(wt1), byrow=TRUE)
#residuals <- fit1$residuals
residuals<-param[wt1,]-t(t(structure(cbind(fit1$coefficients)[fit1$qr$pivot,], names=names(fit1$coefficients))) %*%t(cbind(1,scaled.sumstat[wt1,])))
residuals<-cbind(residuals)
#rrr<<-residuals
the_m<-apply(residuals,FUN=mean,2)
residuals<-sapply(1:numparam,FUN=function(x){residuals[,x]-the_m[x]})
pred <- sapply(1:numparam,FUN=function(x){pred[,x]+the_m[x]})
sigma2<-apply(as.matrix(residuals),FUN=function(x){sum((x)^2 *weights)/sum(weights)},MARGIN=2)
aic<-sum(wt1)*sum(log(sigma2))+2*(nss+1)*numparam
bic<-sum(wt1)*sum(log(sigma2))+log(sum(wt1))*(nss+1)*numparam
if(hcorr == TRUE){
fit2 <- lsfit(scaled.sumstat[wt1,],log(residuals^2),wt=weights)
auxaux<-t(structure(cbind(fit2$coefficients)[fit2$qr$pivot,], names=names(fit2$coefficients))) %*% c(1,target)
pred.sd <- sqrt(exp(auxaux))
pred.sd <- matrix(pred.sd, nrow=sum(wt1), ncol=numparam, byrow=T)
pred.si<-t(t(structure(cbind(fit2$coefficients)[fit2$qr$pivot,], names=names(fit2$coefficients))) %*%t(cbind(1,scaled.sumstat[wt1,])))
pred.si<-sqrt(exp(pred.si))
#residuals2<-log(residuals^2)-t(t(structure(cbind(fit2$coefficients)[fit2$qr$pivot,], names=names(fit2$coefficients))) %*%t(cbind(1,scaled.sumstat[wt1,])))
#fv <- sqrt(exp(log(fit1$residuals^2) - fit2$residuals))
#fv <- sqrt(exp(residuals2))
# adj.values <- pred + (pred.sd*residuals) / fv # correction heteroscedasticity
adj.values <- pred + (pred.sd*residuals) / pred.si
residuals<-(pred.sd*residuals) / pred.si
}
else{
adj.values <- pred + residuals
}
colnames(adj.values) <- colnames(unadj.values)
lambda <- NULL
}
## neural network regression
## ##########################
if(method == "neuralnet"){
linout <- TRUE
## normalise parameters
param.mad <- c()
for(i in 1:numparam){
param.mad[i] <- mad(param[,i][gwt]) # save for renormalisation
param[,i] <- normalise(param[,i],param[,i][gwt])
}
lambda <- sample(lambda, numnet, replace=T)
fv <- array(dim=c(sum(wt1), numparam, numnet))
pred <- matrix(nrow=numparam, ncol=numnet)
for (i in 1:numnet){
fit1 <- nnet(scaled.sumstat[wt1,], param[wt1,], weights = weights, decay = lambda[i],
size = sizenet, trace = trace, linout = linout, maxit = maxit, ...)
cat(i)
fv[,,i] <- fit1$fitted.values
pred[,i] <- predict(fit1, data.frame(rbind(target)))
}
cat("\n")
pred.med <- apply(pred, 1, median)
pred.med <- matrix(pred.med, nrow=sum(wt1), ncol=numparam, byrow=T)
fitted.values <- apply(fv, c(1,2), median)
residuals <- param[wt1,] - fitted.values # median of fitted values par nnets for each accepted point and parameter
if(hcorr == TRUE){
pred2 <- matrix(nrow=numparam, ncol=numnet)
fv2 <- array(dim=c(sum(wt1), numparam, numnet))
for (i in 1:numnet){
fit2 <- nnet(scaled.sumstat[wt1,], log(residuals^2), weights = weights, decay = lambda[i], size = sizenet, trace = trace, linout = linout, ...)
cat(i)
fv2[,,i] <- fit2$fitted.values
pred2[,i] <- predict(fit2, data.frame(rbind(target)))
}
cat("\n")
pred.sd <- sqrt(exp(apply(pred2, 1, median)))
pred.sd <- matrix(pred.sd, nrow=sum(wt1), ncol=numparam, byrow=T)
fv.sd <- sqrt(exp(apply(fv2, c(1,2), median)))
adj.values <- pred.med + (pred.sd*residuals)/fv.sd # correction heteroscedasticity
residuals<-(pred.sd*residuals)/fv.sd
}
else{
adj.values <- pred.med + residuals
}
colnames(adj.values) <- colnames(unadj.values)
## renormalise
for(i in 1:numparam){
## residuals[,i] <- residuals[,i]*param.mad[i] not much sense...
adj.values[,i] <- adj.values[,i]*param.mad[i]
}
} # end of neuralnet
if(method == "ridge"){
#print("la")
## normalise parameters
param.mad <- c()
for(i in 1:numparam){
param.mad[i] <- mad(param[,i][gwt]) # save for renormalisation
param[,i] <- normalise(param[,i],param[,i][gwt])
}
# print("ici")
numnet<-length(lambda)
#lambda <- sample(lambda, numnet, replace=T)
fv <- array(dim=c(sum(wt1), numparam, numnet))
pred <- matrix(nrow=numparam, ncol=numnet)
mataux<-sqrt(diag(weights))
paramaux<-as.matrix(mataux%*%param[wt1,])
#print(dim(paramaux))
scaledaux<-mataux%*%scaled.sumstat[wt1,]
#targetaux<-drop(mataux%*%target
for(parcur in (1:numparam))
{
#cat(parcur)
fit1<-lm.ridge(paramaux[,parcur]~scaledaux,lambda=lambda)
for (i in 1:numnet){
# fit1 <- nnet(scaled.sumstat[wt1,], param[wt1,], weights = weights, decay = lambda[i],
# size = sizenet, trace = trace, linout = linout, maxit = maxit, ...)
# cat(i)
fv[,parcur,i] <- drop(cbind(1,scaled.sumstat[wt1,])%*%(rbind(coef(fit1))[i,]))
pred[parcur,i] <- drop(c(1, target) %*%(rbind(coef(fit1))[i,]))
}
}
# cat("\n")
pred.med <- apply(pred, 1, median)
pred.med <- matrix(pred.med, nrow=sum(wt1), ncol=numparam, byrow=T)
fitted.values <- apply(fv, c(1,2), median)
residuals <- param[wt1,] - fitted.values # median of fitted values par nnets for each accepted point and parameter
if(hcorr == TRUE){
pred2 <- matrix(nrow=numparam, ncol=numnet)
fv2 <- array(dim=c(sum(wt1), numparam, numnet))
for(parcur in (1:numparam))
{
lresidaux<-(mataux%*%(log(residuals[,parcur]^2)))
fit2<-lm.ridge(lresidaux~scaledaux,lambda=lambda)
for (i in 1:numnet){
#cat(i)
fv2[,parcur,i] <- drop(cbind(1,scaled.sumstat[wt1,])%*%(rbind(coef(fit2))[i,]))
pred2[parcur,i] <- drop(c(1, target) %*%(rbind(coef(fit2))[i,]))
}
}
cat("\n")
pred.sd <- sqrt(exp(apply(pred2, 1, median)))
pred.sd <- matrix(pred.sd, nrow=sum(wt1), ncol=numparam, byrow=T)
fv.sd <- sqrt(exp(apply(fv2, c(1,2), median)))
adj.values <- pred.med + (pred.sd*residuals)/fv.sd # correction heteroscedasticity
residuals<- (pred.sd*residuals)/fv.sd
}
else{
adj.values <- pred.med + residuals
}
colnames(adj.values) <- colnames(unadj.values)
## renormalise
for(i in 1:numparam){
## residuals[,i] <- residuals[,i]*param.mad[i] not much sense...
adj.values[,i] <- adj.values[,i]*param.mad[i]
}
} # end of ridge
} # end of else than rejmethod
## back transform parameter values
## ################################
if(numparam == 1){
unadj.values <- matrix(unadj.values, ncol=1)
if(method!="rejection")
{adj.values <- matrix(adj.values, ncol=1)
residuals <- matrix(residuals, ncol=1)}
}
for (i in 1:numparam){
if(transf[i] == "log"){
unadj.values[,i] <- exp(unadj.values[,i])
adj.values[,i] <- exp(adj.values[,i])
}
else if(transf[i] == "logit"){
unadj.values[,i] <- exp(unadj.values[,i])/(1+exp(unadj.values[,i]))
unadj.values[,i] <- unadj.values[,i]*(logit.bounds[i,2]-logit.bounds[i,1])+logit.bounds[i,1]
adj.values[,i] <- exp(adj.values[,i])/(1+exp(adj.values[,i]))
adj.values[,i] <- adj.values[,i]*(logit.bounds[i,2]-logit.bounds[i,1])+logit.bounds[i,1]
}
}
abc.return(transf, logit.bounds, method, call, numparam, nss, paramnames, statnames,
unadj.values, adj.values, ss, weights, residuals, dist, wt1, gwt, lambda, hcorr,aic,bic)
}
abc.return <- function(transf, logit.bounds, method, call, numparam, nss, paramnames, statnames,
unadj.values, adj.values, ss, weights, residuals, dist, wt1, gwt, lambda, hcorr,aic,bic){
if(method == "rejection"){
out <- list(unadj.values=unadj.values, ss=ss, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="rejection", numparam=numparam, numstat=nss, names=list(parameter.names=paramnames, statistics.names=statnames))
}
else if(method == "loclinear"){
out <- list(adj.values=adj.values, unadj.values=unadj.values,
ss=ss, weights=weights, residuals=residuals, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="loclinear", hcorr = hcorr, numparam=numparam, numstat=nss,aic=aic,bic=bic,
names=list(parameter.names=paramnames, statistics.names=statnames))
}
else if(method == "ridge"){
out <- list(adj.values=adj.values, unadj.values=unadj.values,
ss=ss, weights=weights, residuals=residuals, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="ridge", hcorr = hcorr, numparam=numparam, numstat=nss,
names=list(parameter.names=paramnames, statistics.names=statnames))
}
else if(method == "neuralnet"){
out <- list(adj.values=adj.values, unadj.values=unadj.values,
ss=ss, weights=weights, residuals=residuals, dist=dist,
call=call, na.action=gwt, region=wt1, transf=transf, logit.bounds = logit.bounds,
method="neuralnet", hcorr = hcorr, lambda=lambda, numparam=numparam, numstat=nss,
names=list(parameter.names=paramnames, statistics.names=statnames))
}
class(out) <- "abc"
out
}
is.abc <- function(x){
if (inherits(x, "abc")) TRUE
else FALSE
}
print.abc <- function(x, ...){
if (!inherits(x, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- x
cl <- abc.out$call
cat("Call:\n")
dput(cl, control=NULL)
cat("Method:\n")
## rejection
if (abc.out$method == "rejection")
cat("Rejection\n\n")
## loclinear
else if (abc.out$method == "loclinear"){
if(abc.out$hcorr == T){
cat("Local linear regression\n")
cat("with correction for heteroscedasticity\n\n")
}
else cat("Local linear regression\n\n")
}
## ridge
else if (abc.out$method == "ridge"){
if(abc.out$hcorr == T){
cat("Ridge regression\n")
cat("with correction for heteroscedasticity\n\n")
}
else cat("Ridge regression\n\n")
}
## nnet
else if (abc.out$method == "neuralnet"){
if(abc.out$hcorr == T){
cat("Non-linear regression via neural networks\n")
cat("with correction for heteroscedasticity\n\n")
}
else cat("Non-linear regression via neural networks\n\n")
}
cat("Parameters:\n")
cat(abc.out$names$parameter.names, sep=", ")
cat("\n\n")
cat("Statistics:\n")
cat(abc.out$names$statistics.names, sep=", ")
cat("\n\n")
cat("Total number of simulations", length(abc.out$na.action), "\n\n")
cat("Number of accepted simulations: ", dim(abc.out$unadj.values)[1], "\n\n")
invisible(abc.out)
}
summary.abc <- function(object, unadj = FALSE, intvl = .95, print = TRUE, digits = max(3, getOption("digits")-3), ...){
if (!inherits(object, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- object
np <- abc.out$numparam
cl <- abc.out$call
parnames <- abc.out$names[[1]]
npri <- dim(abc.out$unadj.values)[1]
npost <- dim(abc.out$unadj.values)[1]
myprobs <- c((1-intvl)/2, 0.5, 1-(1-intvl)/2)
if(print){
cat("Call: \n")
dput(cl, control=NULL)
}
if(abc.out$method == "rejection" || unadj){
if(print) cat(paste("Data:\n abc.out$unadj.values (",npost," posterior samples)\n\n", sep=""))
res <- matrix(abc.out$unadj.values, ncol=np)
mymin <- apply(res, 2, min)
mymax <- apply(res, 2, max)
mymean <- apply(res, 2, mean)
mymode <- apply(res, 2, getmode, ...)
quants <- apply(res, 2, quantile, probs = myprobs)
sums <- rbind(mymin, quants[1,], quants[2,], mymean, mymode, quants[3,], mymax)
dimnames(sums) <- list(c("Min.:", paste(c(myprobs[1])*100, "% Perc.:", sep=""),
"Median:", "Mean:", "Mode:",
paste(c(myprobs[3])*100, "% Perc.:", sep=""), "Max.:"),
parnames)
}
else{
if(print){
cat(paste("Data:\n abc.out$adj.values (",npost," posterior samples)\n", sep=""))
cat(paste("Weights:\n abc.out$weights\n\n", sep=""))
}
res <- matrix(abc.out$adj.values, ncol=np)
wt <- abc.out$weights
mymin <- apply(res, 2, min)
mymax <- apply(res, 2, max)
mymean <- apply(res, 2, weighted.mean, w = wt)
mymode <- apply(res, 2, getmode, wt/sum(wt), ...)
quants <- apply(res, 2, function(x) rq(x~1, tau = myprobs, weights = wt)$coef)
sums <- rbind(mymin, quants[1,], quants[2,], mymean, mymode, quants[3,], mymax)
dimnames(sums) <- list(c("Min.:", paste("Weighted ", c(myprobs[1])*100, " % Perc.:", sep=""),
"Weighted Median:", "Weighted Mean:", "Weighted Mode:",
paste("Weighted ", c(myprobs[3])*100, " % Perc.:", sep=""), "Max.:"),
parnames)
}
class(sums) <- "table"
if(print){
if(length(digits) == 1)
print(round(sums, digits = digits), quote=FALSE)
else
print(apply(rbind(digits, sums), 2, function(a) round(a[-1], digits = a[1])), quote=FALSE)
invisible(sums)
}
else sums
}
getmode <- function(x, weights = NULL, ...){
## if(missing(bw) | missing(kernel) | missing(window) | missing(n)) warning("density.default() was used to calculate the posterior mode")
d <- density(x, weights = weights)
d$x[which(d$y == max(d$y))][1]
}
hist.abc <- function(x, unadj = FALSE, true = NULL,
file = NULL, postscript = FALSE, onefile = TRUE, ask = !is.null(deviceIsInteractive()),
col.hist = "grey", col.true = "red", caption = NULL, ...){
if (!inherits(x, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- x
np <- abc.out$numparam
parnames <- abc.out$names[[1]]
if(is.null(caption)) mycaption <- as.graphicsAnnot(parnames)
else mycaption <- caption
## checks if true is given
if(!is.null(true)){
if(!is.vector(true)){
stop("Supply true parameter value(s) as a vector.", call.=F)
}
if (length(true) != np){
stop("Number of true values has to be the same as number of parameters in 'x'.", call.=F)
}
cond <- isTRUE(c(match(names(true), parnames), match(names(true), parnames)))
if(cond) stop("Names do not match in 'x' and 'true'.", call.=F)
}
if(abc.out$method == "rejection") res <- abc.out$unadj.values
else if (unadj){
res <- abc.out$unadj.values
}
else res <- abc.out$adj.values
## Devices
## ##########
save.devAskNewPage <- devAskNewPage()
if(!is.null(file)){
file <- substitute(file)
if(!postscript) pdf(file = paste(file, "pdf", sep="."), onefile=onefile)
if(postscript) postscript(file = paste(file, "ps", sep="."), onefile=onefile)
}
else{
if (ask && prod(par("mfcol")) < np) {
devAskNewPage(TRUE)
}
}
myhist <- list()
for(i in 1:np){
myxlim <- range(res[,i], true[i])
myhist[[i]] <- hist(res[,i], xlab = parnames[i], main = paste("Posterior histogram of ", mycaption[i], sep=""),
col = col.hist, xlim = myxlim, ...)
myhist[[i]]$xname <- parnames[i]
if(!is.null(true)){
abline(v = true[i], lwd = 2, col = col.true)
}
}
if(!is.null(file)){
dev.off()
}
else devAskNewPage(save.devAskNewPage)
names(myhist) <- parnames
invisible(myhist)
} # end of hist.abc
##########################################################################################
## plots: for each parameter
## 1. prior density
## 2. posterior + prior density
## 3. distances vs parameters
## 4. histogram of the residuals
##########################################################################################
plot.abc <- function(x, param, subsample = 1000, true = NULL,
file = NULL, postscript = FALSE, onefile = TRUE, ask = !is.null(deviceIsInteractive()), ...){
if (!inherits(x, "abc"))
stop("Use only with objects of class \"abc\".", call.=F)
abc.out <- x
mymethod <- abc.out$method
if(mymethod == "rejection")
stop("Diagnostic plots can be displayed only when method is \"loclinear\", \"neuralnet\" or \"ridge\".", cal.=F)
if(!is.matrix(param) && !is.data.frame(param) && !is.vector(param)) stop("'param' has to be a matrix, data.frame or vector.", call.=F)
if(is.null(dim(param))) param <- matrix(param, ncol=1)
if(is.data.frame(param)) param <- as.matrix(param)
np <- abc.out$numparam
numsim <- length(param)/np
alldist <- log(abc.out$dist)
myregion <- abc.out$region
residuals <- abc.out$residuals
parnames <- abc.out$names$parameter.names
transf <- abc.out$transf
logit.bounds <- abc.out$logit.bounds
##check if param is compatible with x
cond <- isTRUE(c(match(colnames(param), parnames), match(parnames, colnames(param))))
if(cond) stop("'abc.out' and 'param' are not compatible; paramater names are different.", call.=F)
## checks if true is given
if(!is.null(true)){
if(!is.vector(true)){
stop("Supply true parameter value(s) as a vector.", call.=F)
}
if (length(true) != np){
stop("Number of true values has to be the same as number of parameters in 'x'.", call.=F)
}
cond <- isTRUE(c(match(names(true), parnames), match(names(true), parnames)))
if(cond) stop("Names do not match in 'x' and 'true'.", call.=F)
}
rej <- abc.out$unadj.values
res <- abc.out$adj.values
if(np == 1) rej <- matrix(rej, ncol=1)
if(is.vector(param)){
np.orig <- 1
nsim <- length(param)
}
else if(is.matrix(param)){
np.orig <- dim(param)[2]
nsim <- dim(param)[1]
myorder <- match(parnames, colnames(param))
if(isTRUE(myorder-1:np)){
param <- param[, myorder]
warning("'param' is being re-ordered according to 'abc.out'...", call.=F, immediate.=T)
}
}
## check if param has the right dimensions
if(np.orig != np){
stop("The number parameters supplied in \"param\" is different from that in \"x\".", call.=F)
}
## for (i in 1:np){
## if(transf[i] == "log"){
## if(min(param[,i]) <= 0){
## warning("Correcting out of bounds values for \"log\" transformed parameters.", call.=F, immediate.=T)
## x.tmp <- ifelse(param[,i] <= 0,max(param[,i]),param[,i])
## x.tmp.min <- min(x.tmp)
## param[,i] <- ifelse(param[,i] <= 0, x.tmp.min,param[,i])
## }
## param[,i] <- log(param[,i])
## if(min(rej[,i]) <= 0){
## warning("Correcting out of bounds values for \"log\" transformed parameters.", call.=F, immediate.=T)
## x.tmp <- ifelse(rej[,i] <= 0,max(rej[,i]),rej[,i])
## x.tmp.min <- min(x.tmp)
## rej[,i] <- ifelse(rej[,i] <= 0, x.tmp.min,rej[,i])
## }
## rej[,i] <- log(rej[,i])
## if(min(res[,i]) <= 0){
## warning("Correcting out of bounds values for \"log\" transformed parameters.", call.=F, immediate.=T)
## x.tmp <- ifelse(res[,i] <= 0,max(res[,i]),res[,i])
## x.tmp.min <- min(x.tmp)
## res[,i] <- ifelse(res[,i] <= 0, x.tmp.min,res[,i])
## }
## res[,i] <- log(res[,i])
## if(!is.null(true)){
## if(min(true[i]) <= 0){
## x.tmp <- ifelse(true[i] <= 0,max(true[i]),true[i])
## x.tmp.min <- min(x.tmp)
## true[,i] <- ifelse(true[i] <= 0, x.tmp.min,true[i])
## }
## true[i] <- log(true[i])
## }
## }
## else if(transf[i] == "logit"){
## if(min(param[,i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(param[,i] <= logit.bounds[i,1],max(param[,i]),param[,i])
## x.tmp.min <- min(x.tmp)
## param[,i] <- ifelse(param[,i] <= logit.bounds[i,1], x.tmp.min,param[,i])
## }
## if(max(param[,i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(param[,i] >= logit.bounds[i,2],min(param[,i]),param[,i])
## x.tmp.max <- max(x.tmp)
## param[,i] <- ifelse(param[,i] >= logit.bounds[i,2], x.tmp.max,param[,i])
## }
## param[,i] <- (param[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## param[,i] <- log(param[,i]/(1-param[,i]))
## if(min(rej[,i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(rej[,i] <= logit.bounds[i,1],max(rej[,i]),rej[,i])
## x.tmp.min <- min(x.tmp)
## rej[,i] <- ifelse(rej[,i] <= logit.bounds[i,1], x.tmp.min,rej[,i])
## }
## if(max(rej[,i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(rej[,i] >= logit.bounds[i,2],min(rej[,i]),rej[,i])
## x.tmp.max <- max(x.tmp)
## rej[,i] <- ifelse(rej[,i] >= logit.bounds[i,2], x.tmp.max,rej[,i])
## }
## rej[,i] <- (rej[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## rej[,i] <- log(rej[,i]/(1-rej[,i]))
## if(min(res[,i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(res[,i] <= logit.bounds[i,1],max(res[,i]),res[,i])
## x.tmp.min <- min(x.tmp)
## res[,i] <- ifelse(res[,i] <= logit.bounds[i,1], x.tmp.min,res[,i])
## }
## if(max(res[,i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(res[,i] >= logit.bounds[i,2],min(res[,i]),res[,i])
## x.tmp.max <- max(x.tmp)
## res[,i] <- ifelse(res[,i] >= logit.bounds[i,2], x.tmp.max,res[,i])
## }
## res[,i] <- (res[,i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## res[,i] <- log(res[,i]/(1-res[,i]))
## if(!is.null(true)){
## if(min(true[i]) <= logit.bounds[i,1]){
## x.tmp <- ifelse(true[i] <= logit.bounds[i,1],max(true[i]),true[i])
## x.tmp.min <- min(x.tmp)
## true[i] <- ifelse(true[i] <= logit.bounds[i,1], x.tmp.min,true[i])
## }
## if(max(true[i]) >= logit.bounds[i,2]){
## x.tmp <- ifelse(true[i] >= logit.bounds[i,2],min(true[i]),true[i])
## x.tmp.max <- max(x.tmp)
## true[i] <- ifelse(true[i] >= logit.bounds[i,2], x.tmp.max,true[i])
## }
## true[i] <- (true[i]-logit.bounds[i,1])/(logit.bounds[i,2]-logit.bounds[i,1])
## true[i] <- log(true[i]/(1-true[i]))
## }
## }
## } # end of parameter transformations
## devices
## ########
save.devAskNewPage <- devAskNewPage()
if(!is.null(file)){
file <- substitute(file)
if(!postscript) pdf(file = paste(file, "pdf", sep="."), onefile=onefile)
if(postscript) postscript(file = paste(file, "ps", sep="."), onefile=onefile)
}
else{
if (ask && 1 < np) {
devAskNewPage(TRUE)
}
}
## if param is too large, draw a random sample from it the size of
## which can be set by the user
mysample <- sample(1:nsim, subsample)
## plots
## #####
# prior, reg, rej, true
ltys <- c(3, 1, 1, 3)
lwds <- c(1, 2, 1, 2)
cols <- c("black", "red", "black", "black")
unadj <- TRUE
for(i in 1:np){
par(mfcol=c(2,2))
prior.d <- density(param[mysample,i])
post.d <- density(res[,i])
rej.d <- density(rej[,i])
myxlim <- range(c(post.d$x, rej.d$x, true))
myylim <- range(c(post.d$y, rej.d$y, prior.d$y))
par(cex = 1, cex.main = 1.2, cex.lab = 1.1, lwd = 2)
##if(transf[i] == "none")
myxlab <- parnames[i]
##else if(transf[i] == "log") myxlab <- paste("log(", parnames[i], ")", sep="")
##else if(transf[i] == "logit") myxlab <- paste("log(", parnames[i], "/(1-",parnames[i],"))", sep="")
## 1. prior
## ########
plot(prior.d, main = "", col=cols[1], lty = ltys[1], lwd=lwds[1], xlab=myxlab, ...)
title("Prior", sub = paste("N =", prior.d$n, " Bandwidth =", formatC(prior.d$bw)))
## 2. posterior
## ############
plot(post.d, col=cols[2], lty = ltys[2], lwd=lwds[2], xlim = myxlim, ylim = myylim, main = "", xlab=myxlab, ...)
lines(prior.d, col=cols[1], lty = ltys[1], lwd=lwds[1], ...)
lines(rej.d, col = cols[3], lty = ltys[3], lwd=lwds[3], ...)
if(!is.null(true)){
segments(x0 = true[i], x1 = true[i], y0 = myylim[1], y1 = myylim[2], col=cols[4], lty = ltys[4], lwd=lwds[4])
mtext(text="True value", side = 1, line=2, at=true[i])
}
title(paste("Posterior with \"", mymethod, "\"", sep=""), sub=paste("N =", post.d$n, " Bandwidth =", formatC(post.d$bw)))
mtext("\"rejection\" and prior as reference", side=3, line = 0.5)
## 3. distances
## ############
mypch <- 19
myylim <- range(alldist[mysample])#*c(1,1.2)
plot(param[mysample,i], alldist[mysample], col=cols[1],
xlab=myxlab, ylab="log Euclidean distance", ylim=myylim, pch=mypch, ...)
title("Euclidean distances")
mtext(paste("N(All / plotted) = ", numsim, "/", subsample, sep=" "), side = 3, line = .5)
points(param[myregion,i], alldist[myregion], col=cols[2], pch=mypch, ...)
mypar <- par()
if(!is.null(true)){
segments(x0 = true[i], x1 = true[i], y0 = myylim[1], y1 = myylim[2], col=cols[4], lty = ltys[4], lwd=lwds[4])
mtext(text="True value", side = 1, line=2, at=true[i])
}
## 4. residuals
## ############
if(mymethod=="loclinear") mymain <- "Residuals from lsfit()"
if(mymethod=="neuralnet") mymain <- "Residuals from nnet()"
if(mymethod=="ridge") mymain <- "Residuals from lm.ridge()"
qqnorm(residuals, pch=mypch, main=mymain, sub="Normal Q-Q plot", xlab="Theoretical quantiles", ylab="Residuals",...)
qqline(residuals)
} # np
par(mfcol=c(1,1))
if(!is.null(file)){
dev.off()
}
else devAskNewPage(save.devAskNewPage)
invisible()
} # end of plot.abc
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.