Nothing
## unique procedure for all optimization methods
calcProfileLR <- function(varNames=blackbox.getOption("fittedNames"),
pairlist=list(),
cleanResu="") {
INFO <- list(DemographicModel=blackbox.getOption("DemographicModel"),
spec2DProfiles=blackbox.getOption("spec2DProfiles"),
plotOptions=blackbox.getOption("plotOptions"),
#rosglobal=blackbox.getOption("rosglobal"), ## NOTE this cannot be here
# because rosglobal may be modified by this function
fittedparamnbr=blackbox.getOption("fittedparamnbr"),
graphicPars=blackbox.getOption("graphicPars"),
interactiveGraphics=blackbox.getOption("interactiveGraphics")
)
## first construct a list of all desired pairs of variables
if (length(pairlist)==0L) {
if (INFO$fittedparamnbr>2L) { ## if profiles required
if (length(INFO$spec2DProfiles)>0L) { ## if speficic profiles required
for (locit in seq_len(length(INFO$spec2DProfiles))) {
locstring <- INFO$spec2DProfiles[locit]
locstring <- substr(locstring,2,(nchar(locstring)-1L)) ## removes "(" ")"
vars <- strsplit(locstring,",")[[1]]
pairlist <- c(pairlist, list(c(vars[1], vars[2])))
}
} else if ("All2DProfiles" %innc% INFO$plotOptions ) {
for (kk in 1:INFO$fittedparamnbr) for (ll in 1:INFO$fittedparamnbr) {
if (kk<ll) pairlist <- c(pairlist, list(c(varNames[kk], varNames[ll])))
}
} # else pairlist <- c()
## adds 2Nmu, 2Nm if IBD and not already there
if (("IBD" %in% INFO$DemographicModel) &&
! any(as.logical(lapply(pairlist, function(v){all(c("twoNmu", "twoNm")==v)})))) {
pairlist <- c(pairlist, list(c("twoNmu", "twoNm")))
}
## adds 2Nm, g if IBD and not already there ## FR->FR voir a un meilleur controle utilisateur...
if (("IBD" %in% INFO$DemographicModel) && ! any(as.logical(lapply(pairlist, function(v){all(c("twoNm", "g")==v)})))) {
pairlist <- c(pairlist, list(c("twoNm", "g")))
}
## adds 2Nmu, Nb if IBD and not already there
if (("IBD" %in% INFO$DemographicModel) && ! any(as.logical(lapply(pairlist, function(v){all(c("twoNmu", "latt2Ns2")==v)})))) {
pairlist <- c(pairlist, list(c("twoNmu", "latt2Ns2")))
}
## adds 2Nm, condS2 if IBD, condS2 in varNames and not already there
if (("IBD" %in% INFO$DemographicModel) && "condS2" %in% varNames && ! any(as.logical(lapply(pairlist, function(v){all(c("twoNm", "condS2")==v)})))) {
pairlist <- c(pairlist, list(c("twoNm", "condS2")))
}
} else if (INFO$fittedparamnbr==2L) {pairlist <- list(varNames)}
}
if (INFO$fittedparamnbr>2L) {
localmain <- "Profile likelihood ratio"
} else localmain <- "Likelihood ratio";
prevmsglength <- 0
npaires <- length(pairlist);paireIt <- 0
## a check of feasible pairs: same algo as in calcGridRelProfile
checklist <- list()
for(paire in pairlist) {
testvars <- (paire %w/o% varNames) ## should be empty... except in case adhocked below
## ad hoc fix for Nm bound /Nb profile
if ("latt2Ns2" %in% paire && "twoNm" %in% varNames) {testvars <- (testvars %w/o% "latt2Ns2")}
## ad hoc fix for Nb bound/Nm profile
if ("twoNm" %in% paire && "latt2Ns2" %in% varNames) {testvars <- (testvars %w/o% "twoNm")}
## ad hoc fix...
if ("NMratio" %in% paire && "M1" %in% varNames) {testvars <- (testvars %w/o% "NMratio")}
if ("mratio" %in% paire && "M1" %in% varNames) {testvars <- (testvars %w/o% "mratio")}
if ("m1overmu" %in% paire && "M1" %in% varNames) {testvars <- (testvars %w/o% "m1overmu")}
if ("m2overmu" %in% paire && "M2" %in% varNames) {testvars <- (testvars %w/o% "m2overmu")}
if ("Nratio" %in% paire && "twoNmu" %in% varNames) {testvars <- (testvars %w/o% "Nratio")}
if ("Nancratio" %in% paire && "twoNmu" %in% varNames) {testvars <- (testvars %w/o% "Nancratio")}
if ("NactNfounderratio" %in% paire && "twoNmu" %in% varNames) {testvars <- (testvars %w/o% "NactNfounderratio")}
if ("NfounderNancratio" %in% paire && "twoNancmu" %in% varNames) {testvars <- (testvars %w/o% "NfounderNancratio")}
if ("Dgmu" %in% paire && "D" %in% varNames) {testvars <- (testvars %w/o% "Dgmu")}
if ("Tgmu" %in% paire && "T" %in% varNames) {testvars <- (testvars %w/o% "Tgmu")}
if(length(testvars)==0) checklist=c(checklist, list(paire))
}
pairlist <- checklist ## now contains only pairs valid for calcGridRelProfile
if (length(pairlist)>0) {
message.redef(paste("*** Computing ", length(pairlist), " two-dimensional profile confidence plots *** (may be slow)", sep=""))
#write("\n", file=cleanResu)
}
loccex.axis=INFO$graphicPars$cex.axis
## main loop
for (paire in pairlist) {
paireIt <- paireIt+1
stk <- paire[1];stl <- paire[2]
locxstring <- "";locystring <- "";
if (interactive()) {
tmp <- sapply(paire, formatName, format="ASCII")
locmess <- paste("Computing profile grid for (", paste(tmp, collapse=", "), ")")
locmess <- paste(locmess, " (profile #", paireIt, " out of ", npaires, ")", sep="")
cat(locmess, "\n")
}
loclist <- do.call("calcGridRelProfile",list(fixed=paire))
xGrid <- loclist$xGrid;yGrid <- loclist$yGrid;RelLik <- loclist$RelLik;inKrigSpace <- loclist$inKrigSpace
if ( any( ! is.na(RelLik) )) {
rosglobal <- blackbox.getOption("rosglobal") ## latest one
if (min(length(xGrid), length(yGrid))<2) next; ## no variation for one variable
## we want pretty labels so we must determine xat from pretty xlabs values, not xlabs from pretty xat values...
xargs <- maketicks(stk, labels=xGrid, axis=1, maxticks=INFO$graphicPars$xmaxticks)
xat <- xargs$at;xlabs <- xargs$labels;xlegend <- xargs$legend;xph <- xargs$phantomat
if (stk=="latt2Ns2") {
xMLcoord <- rosglobal$latt2Ns2 ## not times GV$Nbfactor since internal coord of the plot are (log) latt2Ns2
} else xMLcoord <- rosglobal$canonVP[stk]
if(islogscale(stk)) {xMLcoord <- log(xMLcoord)}
## same for y
yargs <- maketicks(stl, labels=yGrid, axis=2, maxticks=INFO$graphicPars$ymaxticks)
yat <- yargs$at;ylabs <- yargs$labels;ylegend <- yargs$legend;yph <- yargs$phantomat
if (stl=="latt2Ns2") {
yMLcoord <- rosglobal$latt2Ns2 ## not times GV$Nbfactor since internal coord of the plot are (log) latt2Ns2
} else yMLcoord <- rosglobal$canonVP[stl]
if(islogscale(stl)) {yMLcoord <- log(yMLcoord)}
if (INFO$interactiveGraphics) provideDevice(bbdefaultPars=TRUE) ## one window for _each_ graphic
RelLik[RelLik>1] <- 1
BadExtrapol <- which( ! inKrigSpace, arr.ind=T) ## inKrigSpace is another (no longer global) array created by calcGridRelProfile
cautiousRelLik <- RelLik;cautiousRelLik[BadExtrapol] <- NA
if (nrow(BadExtrapol)>1) {
maxBad <- max(RelLik[BadExtrapol])
usernames <- sapply(paire, formatName, format="ASCII")
if (maxBad==1) {
warningue <- paste("(!) The highest profile likelihood is extrapolated in the ", paste(usernames, collapse=", "), " profile", sep="")
} else if (maxBad>0.05) {
warningue <- paste(" Some high profile likelihoods are extrapolated in the ", paste(usernames, collapse=", "), " profile", sep="")
} else {
warningue <- paste(" Some profile likelihoods are extrapolated in the ", paste(usernames, collapse=", "), " profile", sep="")
}
message.redef(warningue)
write(warningue, file=cleanResu)
}
if ("cautious" %innc% INFO$plotOptions) {
RelLik <- cautiousRelLik
} else if ("Incautious" %innc% INFO$plotOptions) {
cautiousRelLik <- RelLik
} else { ## default
## contour lines show the uncautious and shading shows the cautious. Analogous, though not exactly, to slice plots
}
## each lik is NA either because (x,y) vars are not in krigSpace or that full (x,y,...) point is not in KrigSpace
if (all(is.na(cautiousRelLik))) {
warningue <- paste(" All profile likelihoods are extrapolated in the ", paste(usernames, collapse=", "), " profile", sep="")
cautiousRelLik <- matrix(-1,nrow=nrow(cautiousRelLik),ncol=ncol(cautiousRelLik)) ## FR->FR quick patch for plot fns
}
niveaux <- c(0.001, 0.01, 0.05, seq(1L, 9L, 1L)/10) ## one contour() call => non_overlapping labels
filled.contour(x=xGrid, y=yGrid, cautiousRelLik, xlab=xlegend,
ylab=ylegend, main=localmain, color.palette=spaMM.colors, #topo.colors, ## col and color have different syntaxes
plot.axes={contour(xGrid, yGrid, RelLik, add=T, nlevels=1, levels=niveaux, labcex=INFO$graphicPars$labcex);
if (!is.null(xph)) axis(1, at=xph, labels=rep("", length(xph)), tcl=-0.1);
if (!is.null(yph)) axis(2, at=yph, labels=rep("", length(yph)), tcl=-0.1);
# points(xGrid[BadExtrapol[, 1]], yGrid[BadExtrapol[, 2]], pch=19); ## (NULL OK on points since limits of plot are known)
points(xMLcoord, yMLcoord, pch="+");
axis(1, at=xat, labels=mantissExp(xlabs), cex.axis=loccex.axis);
axis(2, at=yat, labels=mantissExp(ylabs), cex.axis=loccex.axis)})
if (("BWPlots" %innc% INFO$plotOptions) && length(yGrid)>1) {
if (INFO$interactiveGraphics) provideDevice(bbdefaultPars=TRUE) ## one window for _each_ graphic
filled.contour(x=xGrid, y=yGrid, cautiousRelLik, xlab=xlegend,
ylab=ylegend, main=localmain,
col=gray(seq(0.75,1,len = 20)), # gplots::colorpanel(20, "grey75", "white"), #
plot.axes={contour(xGrid, yGrid, RelLik, add=T, nlevels=1, levels=niveaux, labcex=INFO$graphicPars$labcex);
# points(xGrid[BadExtrapol[, 1]], yGrid[BadExtrapol[, 2]], pch=19);
points(xMLcoord, yMLcoord, pch="+");
axis(1, at=xat, labels=mantissExp(xlabs), cex.axis=loccex.axis);
axis(2, at=yat, labels=mantissExp(ylabs), cex.axis=loccex.axis)
}
)
}
} else {
warningue <- paste("(!) No grid point for the ", paste(paire), " profile was within convex hull of Kriged points", sep="")
message.redef(warningue)
write(warningue, file=cleanResu)
warningue <- "Perhaps increase 'gridSteps' so that some grid points fall within the hull." ## dLnLthreshold currently cannot be changed....
message.redef(warningue)
write(warningue, file=cleanResu)
}
} ## end for paire...
plot1DprofFrom2D()
}
## for labels on contours: the default method="flattest"
## does not always show anything, and the others may be ugly
## for contours only without shading: execute the contour(...) above
## only and without the add=T statement.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.