Nothing
calcPredictorOK <- function(FONKgpointls,
minKrigPtNbr=blackbox.getOption("minKrigPtNbr"),
krigmax=NULL,
topmode=FALSE,
rawPlots=TRUE,
cleanResu=""
) {
fittedNames <- blackbox.getOption("fittedNames")
FONKgNames <- blackbox.getOption("FONKgNames")
fittedparamnbr <- blackbox.getOption("fittedparamnbr")
if (is.null(minKrigPtNbr) || minKrigPtNbr==-1) {
minKrigPtNbr <- floor(500**((fittedparamnbr/3)**(1/2))) ## default minimum point nbr for Kriging 36 159 500 1307 3050 6560
} # else if (minKrigPtNbr<=GCVptnbr) {minKrigPtNbr <- GCVptnbr} ## safe, but restricts experimentation
if (fittedparamnbr==1) minKrigPtNbr <- 90 ## default minimum point nbr for Kriging for OnePop 90 = floor(250**((2/3)**(1/2)))
locstring <- paste(" Selection of points for Kriging:", sep="")
message.redef(locstring)
unselectedKgpt <- FONKgpointls ## for figures and second selectFn() call
## selectFn -> findReplicates -> (in case incorrect imput is detected) -> canonizeFromKrig
## -> uses getOption("FONKgLow") to build message
## => FONKgLow must have been set previously by buildFONKgpointls (and will be reset below)
blob <- selectFn(FONKgpointls,
topmode=topmode,
#, #dLnLthreshold default
designRetain=blackbox.getOption("designRetain"),
minPtNbr=minKrigPtNbr, ## il faut bien utiliser la valeur locale !
maxPtNbr=blackbox.getOption("maxKrigPtNbr"),
replicatesInfo=NA)
FONKgpointls <- blob$ptls
FONKgLow <- apply(FONKgpointls[,FONKgNames, drop=FALSE], 2, min)
FONKgUp <- apply(FONKgpointls[, FONKgNames, drop=FALSE], 2, max)
## updating of previous existing vars: (cf buildFONKgpointls or prepareData)
do.call(blackbox.options, list(FONKgpointls=FONKgpointls,FONKgLow=FONKgLow, FONKgUp=FONKgUp))
nuniquerows <- blob$nuniquerows
blackbox.options(RMSy = blob$RMSy)
if (! is.na(blob$warningue)) write(blob$warningue, file=cleanResu) ## writing this at this step prevents writing it at each selectFn call
if(nuniquerows<floor(250**((fittedparamnbr/3)**(1/3)))) {## threshold for warning in cleanResu: 24 90 250 587 1246 2461
write(paste(" Only ", nuniquerows, " unique coordinates selected for Kriging.\n", sep=""), file=cleanResu)
}
#
blackbox.options(hulls=NULL) ## important to make sure that GV$hulls are computed once and only once for a given FONKgpointls
locstring <- paste(" (see min/maxKrigPtNbr setting(Migraine)/argument(R) for controlling this).", sep="")
message.redef(locstring)
if (rawPlots) {
if (length(intersect(blackbox.getOption("DemographicModel"), c("IBD", "OnePop" )))>0) {
rawProfiles(unselectedKgpt=unselectedKgpt, FONKgpointls=FONKgpointls, "both")
} else {
rawProfiles(unselectedKgpt=unselectedKgpt, FONKgpointls=FONKgpointls, "unselected")
rawProfiles(unselectedKgpt=unselectedKgpt, FONKgpointls=FONKgpointls, "selected")
}
}
if ( ! ("extrapolateOutOfHull" %innc% blackbox.getOption("miscOptions"))) {
providefullhull(fittedNames) ## sets .blackbox.data$options$hulls$Kgtotal
}
krigit <- 0
if (is.null(krigmax) || krigmax<0) { ## no explicit value: apply defaults
if (.Platform$r_arch=="i386") {
#do.call(blackbox.options, list(krigmax=2198, kriglength=2000, krigoverlap=500))
krigmax <- 2198
kriglength <- 2000
krigoverlap <- 500
} else { ## x64
#do.call(blackbox.options, list(krigmax=10000, kriglength=6000, krigoverlap=1000))
krigmax <- 10000
kriglength <- 6000
krigoverlap <- 1000
}
} else {
kriglength <- blackbox.getOption("kriglength")
krigoverlap <- blackbox.getOption("krigoverlap")
}
message.redef("\n*** Kriging ***")
ycolname <- blackbox.getOption("ycolname")
flushCSmoothTable() ## although calcPredictorOK() is called only in migraine use so flushing the pointer table is unlikely to be required
if (nrow(FONKgpointls)<=krigmax) {
fitobject <- generatePredictor(,first=1L, last=nrow(FONKgpointls)) ## -> OKrig -> CKrigcoefs
#fitobject$shat.GCV <- NA ## otherwise error on 'summarizing' fit ## (old comment pre 2015)
} else {
fitobject <- list() ## will contain $x, $y, $fitted.values, $blocmin $blocmax, $Kriglist
class(fitobject) <- c("OKriglistplus", class(fitobject))
krigit <- (nrow(FONKgpointls)-krigoverlap)%/%(kriglength-krigoverlap) ## for krigit+1 blocks
fitobject$blocmin <- numeric(krigit)
fitobject$blocmax <- numeric(krigit) ## length=krigit sufficient for krigit+1 blocks
fitlist <- list()
class(fitlist) <- c("OKriglist", class(fitlist))
for (bloc in 0:krigit) {
fit <- generatePredictor(,first=bloc*(kriglength-krigoverlap)+1, last=min(nrow(FONKgpointls), (bloc+1)*kriglength-bloc*krigoverlap))
#fit$shat.GCV <- NA ## otherwise error on 'summarizing' fit
fitobject$blocmin[bloc+1] <- min(fit$x[, 1]);fitobject$blocmax[bloc+1] <- max(fit$x[, 1])
fitlist[[bloc+1]] <- fit
}
#fitlist[[krigit+1]] <- FONKgpointls ## could be convenient for predict(fitobject)
fitobject$Kriglist <- fitlist
fitobject$x <- FONKgpointls[, fittedNames, drop=FALSE]
fitobject$y <- - FONKgpointls[, ycolname]
}
## A single Krig object or a list of such object is stored in the global 'fitobject' above.
## This fitobject should serve as default argument for purefn.
## Some diagnostics...
## next code should be equiv to
## fitobject$fitted.values <- as.matrix(predict(fitobject), ncol=1)
## tmp <- fitobject$fitted.values-fit$y
## when there is a single block
## differences between predicted and observed values
predKriged <- predict(fitobject, testhull=FALSE)
fitobject$fitted.values <- predKriged ## prediction of kriged points
predError <- predKriged + as.vector(FONKgpointls[, ycolname]) #+ car y/-y
blackbox.options(RMSpred=sqrt(mean((predError)^2))) #note replicated coordinates count twice (as they should)
llocalst <- paste("RMS residual error from predictions of smoothed points: ", prettynum(blackbox.getOption("RMSpred")),
"\n (RMS of variance from all pairs, smoothed or not: ", prettynum(blob$RMSy),")")
cat(llocalst, "\n") ## FR->FR for call through migraine, cat() goes in Rout_1 but not on screen. llocalst could deserve to go on screen
blackbox.options(fitobject=fitobject)
invisible(fitobject)
}
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.