Nothing
.sss.valid <- function(x){
el <- list(data="data.frame", coords="data.frame", grid="matrix",
knots="list", W="matrix", contract="logical", regular="logical")
requireNamespace("methods",quietly=TRUE)
exists.slot <- match(names(el),methods::slotNames(x))
if(!any(is.na(exists.slot))){
if(!is.null(x@data)){
if(nrow(x@data)==nrow(x@coords)){
if(!(x@contract && nrow(x@data)==nrow(x@W) && nrow(x@grid)==ncol(x@W)) || !(!x@contract && nrow(x@data)==ncol(x@W) && nrow(x@grid)==nrow(x@W))){
paste("contract logical does not match with W matrix.")
}
} else {
paste("Note! nrow(data) != nrow(coords).")
}
} else {
if(!(x@contract && nrow(x@grid)==ncol(x@W)) || !(!x@contract && nrow(x@grid)==nrow(x@W))){
paste("contract logical does not match with W matrix.")
}
}
if(ncol(x@coords)==2L && ncol(x@grid)==2L && length(x@knots)==2L){
valid = NULL
for(i in 1:length(el)){
requireNamespace("methods",quietly=TRUE)
valid[i] <- any(grepl(paste("^",el[[i]],"$",sep=""),class(methods::slot(x,names(el)[[i]]))))
}
if(all(valid)) TRUE else paste("Slots",paste(names(el)[!valid],sep="",collapse=", "),"not match their class.")
} else {
paste("Number of columns in coords/grid, or length in knots is not 2.")
}
} else {
paste("Slots",paste(names(el)[is.na(exists.slot)],sep="",collapse=", "),"are missing.")
}
}
sss <- setClass("sss",
slots = list(data="data.frame", coords="data.frame", grid="matrix",
knots="list", W="matrix", contract="logical", regular="logical"),
package = "scpm", validity = .sss.valid)
.createGrid <- function(coords, check.regular = TRUE){
if (!is.matrix(coords))
coords <- as.matrix(coords)
if (length(dim(coords)) > 2L || !is.numeric(coords) || is.complex(coords))
stop("cooords must be a numeric, real matrix")
dnx = dimnames(coords)
if(is.null(dnx)) dnx <- vector("list", 2)
o.c = order(coords[,2],coords[,1])
coords = coords[o.c,]
x.grid = as.numeric(.getKnots(coords[,1]))
y.grid = as.numeric(.getKnots(coords[,2]))
xy.grid = list(abcissa = x.grid, ordinate= y.grid)
names(xy.grid) = dnx[[2]]
y.x = matrix(as.numeric(unlist(strsplit(unlist(lapply(y.grid,paste,x.grid,sep=" "))," "))),ncol=2,byrow=TRUE)
grid.s = cbind(y.x[,2],y.x[,1])
colnames(grid.s) = dnx[[2]]
attr.gs = attributes(grid.s)
attr.gs$grid.list = xy.grid
attr.gs$W = .incidenceSpatial(coords,grid.s)
attr.gs$ordering = o.c
attr.gs$regular = nrow(coords)==nrow(grid.s) && all(apply(coords,1,paste,sep="",collapse=":")==apply(grid.s,1,paste,sep="",collapse=":"))
attributes(grid.s) = attr.gs
return(grid.s)
}
.fixLoc <- function(data,smallest=0.00000010,largest=0.00000125,digits=8){
if(class(data)=="sss"){
coords <- data@coords
} else {
coords <- data
}
nc <- colnames(coords)
locsId <- paste(coords[,1],coords[,2],sep=",")
requireNamespace("stats",quietly=TRUE)
locs <- stats::xtabs(~locsId)
locx <- coords[,1]
locy <- coords[,2]
for(i in 1:length(locs)){
isrep <- NULL
for(j in 1:nrow(coords)){
if(locsId[j]==names(locs)[i]){
isrep <- c(isrep,i)
if(length(isrep) < locs[i]){
usign <- runif(2,0,1)
locx[j] <- coords[j,1] + ifelse(usign[1]>0.5,1,-1)*round(runif(1,smallest,largest),digits)
locy[j] <- coords[j,2] + ifelse(usign[2]>0.5,1,-1)*round(runif(1,smallest,largest),digits)
}
}
}
}
locsNewId <- paste(locx,locy,sep=",")
coordsNew <- data.frame(locx,locy)
names(coordsNew) <- nc
if(class(data)=="sss"){
dataOld <- data.frame(data@data,coords)
data <- .asNewsss(dataOld, coords = coordsNew, coords.col = NULL, data.col = NULL, contract = TRUE, ordering = "grid")
} else {
dataOld <- data.frame(coords)
data <- data.frame(dataOld,coordsNew)
}
return(data)
}
"create.sss" <- function(coords, data, ...){
if(missing(coords)) stop("A two-columns matrix or data-frame of coordinates must be specified")
if(missing(data)) data <- data.frame() else data <- as.data.frame(data)
if(is.data.frame(coords))
coords <- as.matrix(coords)
if(ncol(coords)!=2)
stop("coords must be a two-column matrix or data-frame")
if(nrow(data)>0)
if(nrow(data)!=nrow(coords)) stop("coords and data must have the same number of rows")
X <- cbind(coords,data)
coords.col <- 1:2
data.col <- NULL
if(ncol(data)>0)
data.col <- c(1:ncol(data)) + 2
contract <- TRUE
ordering <- "grid"
if(missing(X) || is.null(X)) stop("X must be specified!")
if(!is.data.frame(X) && !is.matrix(X) && !is.list(X)) stop("X must be a data.frame, matrix or list of elements.")
if(is.matrix(X)) X <- as.data.frame(X)
requireNamespace("methods",quietly=TRUE)
if(is.data.frame(X)){
attr.X = attributes(X)
if(is.null(coords)) coords <- attr.X$coords
if(is.null(coords) && (is.null(coords.col) || !is.numeric(coords.col))) stop("If X is a data.frame/matrix without a coords attribute then coords.col numeric vector must be specified.")
if(!is.null(coords.col)) coords <- X[,coords.col]
if(missing(data.col) || is.null(data.col)) data.col <- 1:ncol(X)
data.X <- X
} else {
if((is.null(coords) && is.null(X$coords) && !methods::.hasSlot(X,"coords")) && dim(coords)[2]!=2) stop("If X is a list without a coords elements/slots then the coords n*2 numeric matrix must be specified")
if(is.null(coords) && !is.null(X$coords) && !methods::.hasSlot(X,"coords")) coords <- X$coords
if(is.null(coords) && is.null(X$coords) && !methods::.hasSlot(X,"coords")) coords <- X@coords
if(is.null(X$data) && !!methods::.hasSlot(X,"data")) data.X <- NULL;dnx <- NULL
if(!is.null(X$data) && !!methods::.hasSlot(X,"data")) data.X <- as.data.frame(X$data)
if(is.null(X$data) && !methods::.hasSlot(X,"data")) data.X <- as.data.frame(X@data)
if(!is.null(coords.col)) coords <- as.matrix(coords)[,coords.col]
if(!exists("data.X")) data.X <- NULL
if(!is.null(data.X) && (missing(data.col) || is.null(data.col))) data.col <- 1:ncol(data.X)
if(!exists("data.col")) data.col <- NULL
attr.X = attributes(data.X)
}
if(!is.null(data.col)){
data.X <- data.frame(data.X[,data.col],row.names=attr.X$row.names)
names(data.X) <- attr.X$names[data.col]
}
i.o = apply(cbind(ifelse(rep(is.null(data.X),nrow(coords)),1,data.X),coords),1,paste,sep="",collapse=":")
grid.X = .createGrid(coords, contract)
attr.new = attributes(grid.X)
if(!is.null(data.X) && !is.null(data.col)){
data.X <- data.frame(data.X[attr.new$ordering,],row.names=attr.X$row.names[attr.new$ordering])
names(data.X) <- attr.X$names[data.col]
} else {
data.X <- data.frame(data.X[attr.new$ordering,],row.names=attr.X$row.names[attr.new$ordering])
names(data.X) <- attr.X$names
}
coords.X = coords[attr.new$ordering,]
grid.data = grid.X
if(ordering=="original"){
f.o = apply(cbind(ifelse(rep(is.null(data.X),nrow(coords)),1,data.X),coords.X),1,paste,sep="",collapse=":")
get.back = match(i.o,f.o)
if(!is.null(data.X) && !is.null(data.col)){
data.X <- data.frame(data.X[get.back,],row.names=attr.X$row.names)
names(data.X) <- attr.X$names[data.col]
} else {
data.X <- data.frame(data.X[get.back,],row.names=attr.X$row.names)
names(data.X) <- attr.X$names
}
coords.X = coords.X[get.back,]
if(attr.new$regular){
grid.data = grid.data[get.back,]
}
if(contract) attr.new$W = attr.new$W[get.back,] else attr.new$W = attr.new$W[,get.back]
}
attributes(grid.data) <- attr.new[1:3]
newdata <- methods::new("sss", data = data.frame(data.X), coords = data.frame(coords.X),
grid = as.matrix(grid.data), knots = attr.new$grid.list,
W = attr.new$W, contract = contract, regular = attr.new$regular)
newdata <- .fixLoc(newdata, ...)
return(newdata)
}
is.sss <- function(x){
requireNamespace("methods",quietly=TRUE)
methods::is(x,"sss")
}
.asNewsss <- function(X, coords = NULL, coords.col = NULL, data.col = NULL, contract = TRUE, ordering = c("original","grid"), ...){
if(missing(X) || is.null(X)) stop("X must be specified!")
.defaultArg("ordering","grid")
if(!is.data.frame(X) && !is.matrix(X) && !is.list(X)) stop("X must be a data.frame, matrix or list of elements.")
if(is.matrix(X)) X <- as.data.frame(X)
requireNamespace("methods",quietly=TRUE)
if(is.data.frame(X)){
attr.X = attributes(X)
if(is.null(coords)) coords <- attr.X$coords
if(is.null(coords) && (is.null(coords.col) || !is.numeric(coords.col))) stop("If X is a data.frame/matrix without a coords attribute then coords.col numeric vector must be specified.")
if(!is.null(coords.col)) coords <- X[,coords.col]
if(missing(data.col) || is.null(data.col)) data.col <- 1:ncol(X)
data.X <- X
} else {
if((is.null(coords) && is.null(X$coords) && !methods::.hasSlot(X,"coords")) && dim(coords)[2]!=2) stop("If X is a list without a coords elements/slots then the coords n*2 numeric matrix must be specified")
if(is.null(coords) && !is.null(X$coords) && !methods::.hasSlot(X,"coords")) coords <- X$coords
if(is.null(coords) && is.null(X$coords) && methods::.hasSlot(X,"coords")) coords <- X@coords
if(is.null(X$data) && !methods::.hasSlot(X,"data")) data.X <- NULL;dnx <- NULL
if(!is.null(X$data) && !methods::.hasSlot(X,"data")) data.X <- as.data.frame(X$data)
if(is.null(X$data) && methods::.hasSlot(X,"data")) data.X <- as.data.frame(X@data)
if(!is.null(coords.col)) coords <- as.matrix(coords)[,coords.col]
if(!exists("data.X")) data.X <- NULL
if(!is.null(data.X) && (missing(data.col) || is.null(data.col))) data.col <- 1:ncol(data.X)
if(!exists("data.col")) data.col <- NULL
attr.X = attributes(data.X)
}
if(!is.null(data.col)){
data.X <- data.frame(data.X[,data.col],row.names=attr.X$row.names)
names(data.X) <- attr.X$names[data.col]
}
i.o = apply(cbind(ifelse(rep(is.null(data.X),nrow(coords)),1,data.X),coords),1,paste,sep="",collapse=":")
grid.X = .createGrid(coords, contract)
attr.new = attributes(grid.X)
if(!is.null(data.X) && !is.null(data.col)){
data.X <- data.frame(data.X[attr.new$ordering,],row.names=attr.X$row.names[attr.new$ordering])
names(data.X) <- attr.X$names[data.col]
} else {
data.X <- data.frame(data.X[attr.new$ordering,],row.names=attr.X$row.names[attr.new$ordering])
names(data.X) <- attr.X$names
}
coords.X = coords[attr.new$ordering,]
grid.data = grid.X
if(ordering=="original"){
f.o = apply(cbind(ifelse(rep(is.null(data.X),nrow(coords)),1,data.X),coords.X),1,paste,sep="",collapse=":")
get.back = match(i.o,f.o)
if(!is.null(data.X) && !is.null(data.col)){
data.X <- data.frame(data.X[get.back,],row.names=attr.X$row.names)
names(data.X) <- attr.X$names[data.col]
} else {
data.X <- data.frame(data.X[get.back,],row.names=attr.X$row.names)
names(data.X) <- attr.X$names
}
coords.X = coords.X[get.back,]
if(attr.new$regular){
grid.data = grid.data[get.back,]
}
if(contract) attr.new$W = attr.new$W[get.back,] else attr.new$W = attr.new$W[,get.back]
}
attributes(grid.data) <- attr.new[1:3]
newdata <- methods::new("sss", data = data.frame(data.X), coords = data.frame(coords.X),
grid = as.matrix(grid.data), knots = attr.new$grid.list,
W = attr.new$W, contract = contract, regular = attr.new$regular)
return(newdata)
}
"as.sss" <- function(X, coords, coords.col, data.col, ...){
if(missing(coords)){coords <- NULL}
if(missing(coords.col)){coords.col <- NULL}
if(missing(data.col)){data.col <- NULL}
newdata <- .asNewsss(X = X, coords = coords, coords.col = coords.col, data.col = data.col, contract = TRUE, ordering = "grid")
newdata <- .fixLoc(newdata, ...)
return(newdata)
}
"sss2df" <- function(x){
if(class(x)=="sss"){
datum <- data.frame(x@coords,x@data)
attr2df <- attributes(x)
attr2df$class <- "data.frame"
attr2df$contract <- NULL
attr2df$type <- NULL
attr(datum,"data") <- attr2df$data
attr(datum,"coords") <- attr2df$coords
attr(datum,"grid") <- attr2df$grid
attr(datum,"knots") <- attr2df$knots
attr(datum,"W") <- attr2df$W
attr(datum,"regular") <- attr2df$regular
return(datum)
} else {
stop("Not an element from class sss")
}
}
linear <- function(formula, data = NULL, contrasts = NULL, intercept = FALSE){
requireNamespace("stats", quietly = TRUE)
intercept <- any(grepl("(^1|+1)",formula))
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "contrasts"), names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
o = match(c("formula"),names(mf))
names(mf)[o] = "object"
mf[[1L]] <- quote(stats::model.matrix)
mf <- eval(mf, parent.frame())
if(!intercept){
if(!is.na(match("(Intercept)",colnames(mf)))){
nn <- colnames(mf)
mf = as.matrix(mf[,-match("(Intercept)",nn)])
colnames(mf) <- nn[-match("(Intercept)",nn)]
}
}
return(mf)
}
cp <- function(x, psi, data = NULL, groups = NULL, contrasts = NULL, only.UV = FALSE){
requireNamespace("stats", quietly = TRUE)
mf <- match.call(expand.dots = FALSE)
m <- match(c("x", "data"), names(mf), 0L)
c0 <- mf[c(1L, m)]
o = match(c("x"),names(c0))
names(c0)[o] = "object"
c0[[1L]] <- quote(stats::model.matrix)
cv <- eval(c0, parent.frame())
if(!is.na(match("(Intercept)",colnames(cv)))){
cv = cv[,-1]
}
if(is.null(groups)){
uv = .matrixUV(cv, psi = psi)
} else {
m <- match(c("groups", "data","contrasts"), names(mf), 0L)
c1 <- mf[c(1L, m)]
if(class(groups)=="formula"){
atg <- attributes(stats::terms(groups))
if(atg$intercept==1){
c1[[2L]] = stats::formula(paste("~ -1 +",paste(atg$term.labels,sep="",collapse="+")))
}
o = match(c("groups"),names(c1))
names(c1)[o] = "object"
c1[[1L]] <- quote(stats::model.matrix)
cc <- eval(c1, parent.frame())
cc = apply(cc,1,paste,sep="")
} else {
cc = groups
}
cn = .getKnots(cc)
cc = .incidenceMatrix(cc)
colnames(cc) = cn
if(ncol(cc)!=length(psi)) stop("Number of elements in psi must be equal to the number of groups in factor.")
cv = apply(cc,2,function(x){return(x*cv)})
uv = NULL
for(i in 1:length(psi)){
uv = cbind(uv,.matrixUV(cv[,i], psi = psi[i]))
}
}
list(X = if(!only.UV){cv}, C = if(!is.null(groups) & !only.UV){cc}, UV = uv, psi = psi, call = mf)
}
s2D = function(data = NULL, penalty = c("none","cs","ps","tps"), is.X = c("none","tensor","tps"), intercept = TRUE,
ps.order = 2, Ztr = c("RVA","RV","VA","V"), Amethod = c("eigen","svd"),
Adiag = TRUE, aniso.angle = 0, aniso.ratio = 1, env = .GlobalEnv, ...){
.defaultArg("penalty","none")
.defaultArg("Ztr","RV")
.defaultArg("Amethod","eigen")
if(missing(is.X)) is.X <- NULL
if(penalty=="none"){
if(is.null(is.X))
stop("If penalty=\"none\" then is.X must be specified and different than \"none\"!")
} else if(penalty=="cs" | penalty=="ps"){
if(is.null(is.X)){
is.X = "tensor"
} else if(is.X!="tensor"){
is.X = "tensor"
message("Changing is.X. If penalty=\"cs\" or \"ps\" then is.X must be \"tensor\"!.")
}
} else if(penalty=="tps"){
if(is.null(is.X)){
is.X = "tps"
} else if(is.X!="tps"){
is.X = "tps"
message("Changing is.X. If penalty=\"tps\" then is.X must be \"tps\"!.")
}
}
if(missing(data)) data <- NULL
if(is.null(data)){
ol = ls(envir = env)
ol.id = match(c("data"),ol)
if(all(is.na(ol.id))){
stop("Input a sss (data) object.")
} else {
if(!is.na(match("data",ol))){
data = get("data", envir = env)
}
}
} else if(!is.sss(data)){
stop("Input a sss (data) object.")
}
if(penalty=="none") smooth <- FALSE else smooth <- TRUE
if(smooth){
if(penalty=="ps" | penalty=="cs"){
c.Q = .tensorBase(obj = data, grid.list = attributes(data@grid)$grid.list, type = penalty, diff.order = ps.order)
} else if(penalty=="tps"){
c.Q = .tpsBase(data@grid, ...)
}
} else {
c.Q <- NULL
}
if(is.X=="none"){
c.X = NULL
} else if(is.X=="tensor" | is.X=="tps"){
c.xym = data@coords
c.X = as.matrix(cbind(if(intercept){1},as.matrix(c.xym)))
colnames(c.X) = c(if(intercept){"Intercept"},colnames(as.matrix(c.xym)))
if(is.X=="tensor"){
aux.cX <- colnames(c.X)
c.X = as.matrix(cbind(c.X,apply(c.xym,1,prod)))
colnames(c.X) = c(aux.cX,paste(colnames(c.xym),sep="",collapse="."))
}
}
c.h = .h2D(as.matrix(data@grid), angle = aniso.angle, ratio = aniso.ratio)$h
c.use = list(penalty = penalty, is.X = is.X, intercept = intercept, ps.order = ps.order, Ztr = Ztr, Amethod = Amethod, Adiag = ifelse(penalty=="ps", FALSE, Adiag), aniso.angle = aniso.angle, aniso.ratio = aniso.ratio)
list(X = c.X, h = c.h, Q = c.Q, use = c.use)
}
s2D = function(data = NULL, penalty = c("none","cs","ps","tps"), is.X = c("none","tensor","tps"), intercept = TRUE,
ps.order = 2, aniso.angle = 0, aniso.ratio = 1, env = .GlobalEnv, ...){
.defaultArg("penalty","none")
Ztr <- "V"
Amethod <- "eigen"
Adiag = TRUE
if(missing(is.X)) is.X <- NULL
if(penalty=="none"){
if(is.null(is.X))
stop("If penalty=\"none\" then is.X must be specified and different than \"none\"!")
} else if(penalty=="cs" | penalty=="ps"){
if(is.null(is.X)){
is.X = "tensor"
} else if(is.X!="tensor"){
is.X = "tensor"
message("Changing is.X. If penalty=\"cs\" or \"ps\" then is.X must be \"tensor\"!.")
}
} else if(penalty=="tps"){
if(is.null(is.X)){
is.X = "tps"
} else if(is.X!="tps"){
is.X = "tps"
message("Changing is.X. If penalty=\"tps\" then is.X must be \"tps\"!.")
}
}
if(missing(data)) data <- NULL
if(is.null(data)){
ol = ls(envir = env)
ol.id = match(c("data"),ol)
if(all(is.na(ol.id))){
stop("Input a sss (data) object.")
} else {
if(!is.na(match("data",ol))){
data = get("data", envir = env)
}
}
} else if(!is.sss(data)){
stop("Input a sss (data) object.")
}
if(penalty=="none") smooth <- FALSE else smooth <- TRUE
l.Z <- NULL
if(smooth){
if(penalty=="ps" | penalty=="cs"){
c.Q = .tensorBase(obj = data, grid.list = attributes(data@grid)$grid.list, type = penalty, diff.order = ps.order)
l.Z <- list(c.Q$Qx$V%*%solve(t(c.Q$Qx$V)%*%c.Q$Qx$V), c.Q$Qy$V%*%solve(t(c.Q$Qy$V)%*%c.Q$Qy$V))
names(l.Z) <- names(data@coords)
} else if(penalty=="tps"){
c.Q = .tpsBase(data@grid, ...)
}
} else {
c.Q <- NULL
}
l.X <- NULL
if(is.X=="none"){
c.X = NULL
} else if(is.X=="tensor" | is.X=="tps"){
c.xym = data@coords
c.X = as.matrix(cbind(if(intercept){1},as.matrix(c.xym)))
colnames(c.X) = c(if(intercept){"Intercept"},colnames(as.matrix(c.xym)))
if(is.X=="tensor"){
aux.cX <- colnames(c.X)
c.X = as.matrix(cbind(c.X,apply(c.xym,1,prod)))
colnames(c.X) = c(aux.cX,paste(colnames(c.xym),sep="",collapse="."))
}
l.X <- list(cbind(if(intercept){1},.getKnots(data@coords[,1])), cbind(if(intercept){1},.getKnots(data@coords[,2])))
names(l.X) <- names(data@coords)
}
l.F <- NULL
if(smooth){
if(penalty!="tps"){
l.F <- list(X = l.X[[2]]%x%l.X[[1]], Z = cbind(l.X[[2]]%x%l.Z[[1]], l.Z[[2]]%x%l.X[[1]], l.Z[[2]]%x%l.Z[[1]]))
}
}
c.h = .h2D(as.matrix(data@grid), angle = aniso.angle, ratio = aniso.ratio)$h
c.use = list(penalty = penalty, is.X = is.X, intercept = intercept, ps.order = ps.order, Ztr = Ztr, Amethod = Amethod, Adiag = ifelse(penalty=="ps", FALSE, Adiag), aniso.angle = aniso.angle, aniso.ratio = aniso.ratio)
list(F = l.F, X = c.X, h = c.h, Q = c.Q, use = c.use)
}
.scp <- function(formula, data,
initial = NULL, contrasts = NULL,
pars = list(fix.nugget = FALSE, fix.kappa = FALSE, is.log.rhos = TRUE,
model = "exponential", nugget.tol = 1.0e-15),
method = list(use.reml = FALSE, use.profile = TRUE, linearize = "eigen"),
control.optim = list(), control.gs = list(), control.ch = list(), print.pars = TRUE){
if(missing(pars) | is.null(pars)){
pars = list(fix.nugget = NULL, fix.kappa = NULL, is.log.rhos = TRUE,
model = "matern", nugget.tol = 1.0e-15)
}
if(missing(method) | is.null(method)){
method = list(use.reml = FALSE, use.profile = TRUE, linearize = "eigen")
}
if(!is.null(pars$fix.nugget)){
if(is.logical(pars$fix.nugget)){
if(pars$fix.nugget){
warning("Not value set for fix.nugget!!\n\nTo fix the nugget to a value different than zero set fix.nugget=value.\nTo estimate the nugget set fix.nugget=FALSE or NULL.")
pars$fix.nugget = 0
} else {
pars$fix.nugget = NULL
}
}
}
if(!is.null(pars$fix.kappa)){
if(is.logical(pars$fix.kappa)){
if(all(pars$fix.kappa)){
warning("Not value set for fix.kappa!!\n\nTo fix kappa to a value different than 0.5 set fix.kappa=value.\nTo estimate kappa set fix.kappa=FALSE or NULL.")
pars$fix.kappa = rep(0.5,length(pars$fix.kappa))
} else {
pars$fix.kappa = NULL
}
}
}
call.v = match.call()
requireNamespace("stats", quietly = TRUE)
mf = stats::terms(formula, specials=c("linear","cp","s1D","s2D"))
af = attributes(mf)
cr = rownames(af$factors)
wl = grepl("linear\\(",cr)
if(sum(wl)==1){
cl = as.pairlist(parse(text=cr[af$specials$linear]))[[1]]
if(is.null(cl$data)) cl$data <- data@data
if(is.null(cl$contrasts)) cl$contrasts <- contrasts
el = eval(cl, parent.frame())
} else el = NULL
wcp = grepl("cp\\(",cr)
if(sum(wcp)>0){
ec = NULL
for(i in 1:sum(wcp)){
cc = as.pairlist(parse(text=cr[af$specials$cp][i]))[[1]]
if(is.null(cc$data)) cc$data <- data@data
if(is.null(cc$contrasts)) cc$contrasts <- contrasts
cc$only.UV <- FALSE
ec[[i]] = eval(cc, parent.frame())
}
} else ec = NULL
ws1D = grepl("s1D\\(",cr)
if(sum(ws1D)>0){
eg = NULL
for(i in 1:sum(ws1D)){
cg = as.pairlist(parse(text=cr[af$specials$s1D]))[[1]]
if(is.null(cg$data)) cg$data <- data@data
eg[[i]] = eval(cg, parent.frame())
}
} else eg = NULL
ws2D = grepl("s2D\\(",cr)
if(sum(ws2D)==1){
cs = as.pairlist(parse(text=cr[af$specials$s2D]))[[1]]
if(!is.null(el)){
if(any(grepl("\\(Intercept\\)",colnames(el)))){
if(ncol(el)>1){
nel <- colnames(el)
el <- as.matrix(el[,-match("(Intercept)",nel)])
colnames(el) <- nel[-match("(Intercept)",nel)]
cs$intercept <- TRUE
} else {
el <- NULL
}
}
}
if(is.null(cs$data)) cs$data <- data
if(is.null(cs$aniso.angle)){ if(!is.null(pars$aniso.angle)){ cs$aniso.angle <- pars$aniso.angle } else { cs$aniso.angle <- 0 }} else { if(is.null(pars$aniso.angle)){ pars$aniso.angle <- cs$aniso.angle } }
if(is.null(cs$aniso.ratio)){ if(!is.null(pars$aniso.ratio)){ cs$aniso.ratio <- pars$aniso.ratio } else { cs$aniso.ratio <- 1 }} else { if(is.null(pars$aniso.ratio)){ pars$aniso.ratio <- cs$aniso.ratio } }
es = eval(cs, parent.frame())
} else {
es = list()
es$use = list(penalty = "none", is.X = "none")
if(is.null(pars$aniso.angle)) pars$aniso.angle <- 0
if(is.null(pars$aniso.ratio)) pars$aniso.ratio <- 1
if(is.null(pars$inverse)) pars$inverse <- "solve"
if(is.null(pars$tol)) pars$tol <- .Machine$double.neg.eps
es$use$inverse <- pars$inverse
es$use$tol <- pars$tol
es$h = .h2D(as.matrix(data@grid), angle = pars$aniso.angle, ratio = pars$aniso.ratio)$h
es$F = list(X = NULL, Z = NULL)
es$X = NULL
es$Q = NULL
}
zres = stats::formula(paste(cr[1],"~ 1"))
cz = match.call(expand.dots = FALSE)
m <- match(c("formula", "data"), names(cz), 0L)
cz <- cz[c(1L, m)]
cz$data <- data@data
cz$formula <- zres
cz[[1L]] <- quote(stats::model.frame)
zV <- eval(cz, parent.frame())
zV <- stats::model.response(zV)
XM = cbind(el)
VM <- NULL
ZM <- NULL
if(!is.null(es$F)){
if(!is.null(es$F$Z)){
ZM = cbind(es$F$Z)
VM <- .generateV(ZM, itol = pars$tol)
}
}
cpc = NULL
cppsi = NULL
UVel = NULL
if(length(ec)>0){
for(i in 1:length(ec)){
ec.aux = ifelse(is.null(XM),0,ncol(XM))
XM = cbind(XM,ec[[i]]$UV)
cpc[[i]] = ec[[i]]$call
cppsi[[i]] = ec[[i]]$psi
UVel[[i]] = (ec.aux+1):ncol(XM)
}
ch.l = list(cp.call = cpc, psi0 = cppsi, UVcols = UVel)
} else {
ch.l = NULL
}
for(i in 1:length(eg)){
XM = cbind(XM,eg[[i]]$X)
}
if(!is.null(XM)){
XM = cbind(XM,es$X)
} else {
XM = cbind(es$X)
}
message("Starting computation")
fit = .scpInitialize(zV = zV, XM = XM, ZM = ZM, VM = VM, hM = es$h, Qel = es$Q, W = data@W, ch = ch.l,
use = es$use, pars = pars, initial = initial,
method = method, control.optim = control.optim,
control.gs = control.gs, control.ch = control.ch, print.pars = print.pars)
message("Ending computation")
message("Organising output")
fit$penalty$matrices = es$Q
out <- sssFit()
out@data <- data
out@zV <- zV
out@XL <- ec
out@XC <- ec
out@XF <- eg
out@XS <- es
out@fit <- fit
out@call <- call.v
out
}
scp <- function(formula, data, initial = NULL, contrasts = NULL,
model = "exponential", fix.nugget = FALSE, fix.kappa = FALSE,
nugget.tol = 1.0e-15, angle = 0, ratio = 1,
use.reml = FALSE, use.profile = TRUE, chMaxiter = 20, control = list()){
mc <- match.call()
if(missing(model)){ model <- "exponential" } else { if(is.null(model)){ model <- "exponential" } }
pars <- list(fix.nugget = fix.nugget, fix.kappa = fix.kappa, is.log.rhos = TRUE, model = model, nugget.tol = nugget.tol, aniso.angle = angle, aniso.ratio = ratio)
if(!is.null(control$tol))
pars$tol <- control$tol
method <- list(use.reml = use.reml, use.profile = use.profile)
control.gs <- control[match(c("is.zero", "is.inf","size.each.grid","size.sample"),names(control))]
control.ch <- list(maxiter = chMaxiter, maxskip = 10)
control.optim <- control[match(c("trace", "fnscale","parscale","ndeps","maxit","abstol","reltol","alpha","beta","gamma","REPORT","type","lmm","factr","pgtol","temp","tmax"),names(control))]
out <- .scp(formula, data, initial, contrasts, pars = pars, method = method, control.optim = control.optim, control.gs = control.gs, control.ch = control.ch, print.pars = FALSE)
out@call <- mc
out
}
.scpLinearApprox <- function(obj,c,proj.method=c("Schur","proj","H-proj","ginv"),proj.to=c("sxz","sz")){
.defaultArg("proj.method","H-proj")
.defaultArg("proj.to","sz")
X0 <- obj$fit$X
R0 <- obj$fit$R
Sigma <- obj$fit$sigma2s%x%R0
H0 <- solve(R0)
P0 <- X0%*%solve(t(X0)%*%H0%*%X0)%*%t(X0)%*%H0
In <- diag(1,nrow(P0))
Vdim <- ifelse(obj$XS$use$is.X=="tensor",4,ifelse(obj$XS$use$is.X=="tps",3,0))
CM <- X0[,-c((ncol(X0)-(Vdim-1)):ncol(X0))]
Xs <- X0[,(ncol(X0)-Vdim+1):ncol(X0)]
Zs <- obj$fit$Zs
pz <- Zs%*%solve(t(Zs)%*%H0%*%Zs + sum(obj$fit$alpha)%x%obj$fit$penalty$decompose$Ai)%*%t(Zs)%*%H0
ipz <- In - pz
px <- Xs%*%solve(t(Xs)%*%H0%*%ipz%*%Xs)%*%t(Xs)%*%H0%*%ipz
ipx <- In - px
sxz <- px + pz%*%ipx
requireNamespace("Matrix",quietly=TRUE)
if(proj.to=="sz"){
ss <- Matrix::Schur(pz)
oo <- order(ss$EValues,decreasing=TRUE)
Qs <- ss$Q[,oo]
Qsc <- Qs[,1:c]
if(proj.method=="Schur")
SS <- Qsc%*%t(Qsc)
if(proj.method=="proj")
SS <- Qsc%*%solve(t(Qsc)%*%Qsc)%*%t(Qsc)
if(proj.method=="H-proj")
SS <- Qsc%*%solve(t(Qsc)%*%H0%*%Qsc)%*%t(Qsc)%*%H0
if(proj.method=="ginv")
SS <- Qsc%*%.whichInverse("ginv",Qsc)
SZ <- SS
ISZ <- In - SZ
SX <- Xs%*%solve(t(Xs)%*%H0%*%ISZ%*%Xs)%*%t(Xs)%*%H0%*%ISZ
ISX <- In - SX
SS <- SX + SZ%*%ISX
} else {
SX <- NULL
SZ <- NULL
ss <- Matrix::Schur(sxz)
oo <- order(ss$EValues,decreasing=TRUE)
Qs <- ss$Q[,oo]
Qsc <- Qs[,1:c]
if(proj.method=="Schur")
SS <- Qsc%*%t(Qsc)
if(proj.method=="proj")
SS <- Qsc%*%solve(t(Qsc)%*%Qsc)%*%t(Qsc)
if(proj.method=="H-proj")
SS <- Qsc%*%solve(t(Qsc)%*%H0%*%Qsc)%*%t(Qsc)%*%H0
if(proj.method=="ginv")
SS <- Qsc%*%.whichInverse("ginv",Qsc)
}
if(ncol(CM)>0){
PS <- CM%*%solve(t(CM)%*%H0%*%(In-SS)%*%CM)%*%t(CM)%*%H0%*%(In-SS)
MS <- PS + SS%*%(In-PS)
} else {
PS <- NULL
MS <- SS
}
B0 <- t(In-P0)%*%solve(Sigma)%*%(In-P0)
SS0 <- t(obj$zV)%*%B0%*%obj$zV
df0 <- sum(diag(B0%*%Sigma))
B1 <- t(In-MS)%*%solve(Sigma)%*%(In-MS)
SS1 <- t(obj$zV)%*%B1%*%obj$zV
df1 <- sum(diag(B1%*%Sigma))
criteria <- (df0 >= df1)
null <- ifelse(criteria,0,1)
if(criteria){
B01 <- B0-B1
} else {
B01 <- B1-B0
}
SS01 <- t(obj$zV)%*%B01%*%obj$zV
df01 <- sum(diag(B01%*%Sigma))
Ft <- ifelse(criteria,(df1/df01)*(SS01/SS1),(df0/df01)*(SS01/SS0))
pv <- 1 - pf(Ft,df01,ifelse(criteria,df1,df0))
requireNamespace("mvtnorm",quietly=TRUE)
loglik0 <- mvtnorm::dmvnorm(as.vector(obj$zV-P0%*%obj$zV),mean=rep(0,length(obj$zV)),sigma=Sigma,log=TRUE)
aic0 <- -2 * loglik0 - 2 * sum(diag(P0))
bic0 <- -2 * loglik0 - sum(diag(P0)) * log(length(obj$zV))
logliks <- mvtnorm::dmvnorm(as.vector(obj$zV-MS%*%obj$zV),mean=rep(0,length(obj$zV)),sigma=Sigma,log=TRUE)
aics <- -2*logliks - 2*sum(diag(SS))
bics <- -2*logliks - sum(diag(SS))*log(length(obj$zV))
o <- list(c=c,null=null,"tr(M0)"=sum(diag(P0)),"tr(S(phi))"=sum(diag(obj$fit$S)),"tr(S(c))"=sum(diag(SS)),"tr(M(phi))"=sum(diag(obj$fit$M)),"tr(M(c))"=sum(diag(MS)),F=round(Ft,5),p=round(pv,6),loglik0=loglik0,aic0=aic0,bic0=bic0,loglik=logliks,aic=aics,bic=bics,df0=df0,df1=round(df1,5),df01=round(df01,5),SS0=round(SS0,5),SS1=round(SS1,5),SS01=round(SS01,5))
attr(o,"SX") <- SX
attr(o,"SZ") <- SZ
attr(o,"PX") <- px
attr(o,"PZ") <- pz
attr(o,"PS") <- PS
attr(o,"A") <- Qsc
attr(o,"Sigma") <- Sigma
attr(o,"M0") <- P0
attr(o,"P1") <- SS
attr(o,"M1") <- MS
attr(o,"B0") <- B0
attr(o,"B1") <- B1
attr(o,"B01") <- B01
attr(o,"H0") <- ifelse(null==0,"Linear","Spline")
attr(o,"H1") <- ifelse(null==0,"Spline","Linear")
return(o)
}
.scpLinearApprox <- function (object, c, proj.method = c("Schur", "proj", "H-proj", "ginv"), proj.to = c("sxz", "sz"), sigma = NULL, R = NULL, tol = .Machine$double.neg.eps*1.0e-10){
.defaultArg("proj.method", "H-proj")
.defaultArg("proj.to", "sz")
if(is.null(object@zV)){
aux <- object
object <- sssFit()
object@zV <- aux$y
object@fit <- list()
object@fit$X <- aux$Xss
object@fit$R <- aux$R
object@fit$sigma2s <- aux$sigma2
object@fit$W <- aux$W
object@fit$Zs <- aux$W%*%aux$Z
object@fit$alpha <- aux$alpha
object@fit$psi <- aux$psi
object@fit$M <- aux$M
object@fit$S <- aux$Sa
object@fit$penalty <- list()
object@fit$penalty$decompose <- list()
object@fit$penalty$decompose$Ai <- solve(.csBase(.getKnots(object@fit$X[,ncol(object@fit$X)]))$A,tol = tol)
object@XS <- list()
object@XS$use <- list()
object@XS$use$is.X <- "uni"
}
if(!is.null(sigma) & !is.null(R)){
object@fit$sigma2s <- sigma**2
object@fit$R <- R
}
X0 <- object@fit$X
R0 <- object@fit$R
Sigma <- object@fit$sigma2s %x% R0
H0 <- solve(R0,tol = tol)
P0 <- X0 %*% solve(t(X0) %*% H0 %*% X0,tol = tol) %*% t(X0) %*% H0
In <- diag(1, nrow(P0))
Vdim <- ifelse(object@XS$use$is.X == "tensor", 4, ifelse(object@XS$use$is.X == "tps", 3, ifelse(object@XS$use$is.X == "uni", 2, 0)))
CM <- cbind(X0[, -c((ncol(X0) - (Vdim - 1)):ncol(X0))])
if(object@XS$use$is.X != "uni"){
Xs <- X0[, (ncol(X0) - Vdim + 1):ncol(X0)]
} else {
Xs <- X0[, (ncol(X0) - Vdim + 1):ncol(X0)]
}
Zs <- object@fit$Zs
pz <- Zs %*% solve(t(Zs) %*% H0 %*% Zs + sum(object@fit$alpha) %x% object@fit$penalty$decompose$Ai,tol = tol) %*% t(Zs) %*% H0
ipz <- In - pz
px <- Xs %*% solve(t(Xs) %*% H0 %*% ipz %*% Xs,tol = tol) %*% t(Xs) %*% H0 %*% ipz
ipx <- In - px
sxz <- px + pz %*% ipx
requireNamespace("Matrix", quietly = TRUE)
if (proj.to == "sz") {
ss <- Matrix::Schur(pz)
oo <- order(ss$EValues, decreasing = TRUE)
Qs <- ss$Q[, oo]
Qsc <- Qs[, 1:c]
if (proj.method == "Schur")
SS <- Qsc %*% t(Qsc)
if (proj.method == "proj")
SS <- Qsc %*% solve(t(Qsc) %*% Qsc,tol = tol) %*% t(Qsc)
if (proj.method == "H-proj")
SS <- Qsc %*% solve(t(Qsc) %*% H0 %*% Qsc,tol = tol) %*% t(Qsc) %*% H0
if (proj.method == "ginv")
SS <- Qsc %*% .whichInverse("ginv", Qsc)
SZ <- SS
ISZ <- In - SZ
SX <- Xs %*% solve(t(Xs) %*% H0 %*% ISZ %*% Xs,tol = tol) %*% t(Xs) %*% H0 %*% ISZ
ISX <- In - SX
SS <- SX + SZ %*% ISX
}
else {
SX <- NULL
SZ <- NULL
ss <- Matrix::Schur(sxz)
oo <- order(ss$EValues, decreasing = TRUE)
Qs <- ss$Q[, oo]
Qsc <- Qs[, 1:c]
if (proj.method == "Schur")
SS <- Qsc %*% t(Qsc)
if (proj.method == "proj")
SS <- Qsc %*% solve(t(Qsc) %*% Qsc,tol = tol) %*% t(Qsc)
if (proj.method == "H-proj")
SS <- Qsc %*% solve(t(Qsc) %*% H0 %*% Qsc,tol = tol) %*% t(Qsc) %*% H0
if (proj.method == "ginv")
SS <- Qsc %*% .whichInverse("ginv", Qsc)
}
if (ncol(CM) > 0) {
PS <- CM %*% solve(t(CM) %*% H0 %*% (In - SS) %*% CM,tol = tol) %*% t(CM) %*% H0 %*% (In - SS)
MS <- PS + SS %*% (In - PS)
}
else {
PS <- NULL
MS <- SS
}
B0 <- t(In - P0) %*% solve(Sigma,tol = tol) %*% (In - P0)
SS0 <- t(object@zV) %*% B0 %*% object@zV
df0 <- sum(diag(B0 %*% Sigma))
B1 <- t(In - MS) %*% solve(Sigma,tol = tol) %*% (In - MS)
SS1 <- t(object@zV) %*% B1 %*% object@zV
df1 <- sum(diag(B1 %*% Sigma))
criteria <- (df0 >= df1)
null <- ifelse(criteria, 0, 1)
if (criteria) {
B01 <- B0 - B1
}
else {
B01 <- B1 - B0
}
SS01 <- t(object@zV) %*% B01 %*% object@zV
df01 <- sum(diag(B01 %*% Sigma))
Ft <- ifelse(criteria, (df1/df01) * (SS01/SS1), (df0/df01) * (SS01/SS0))
pv <- 1 - pf(Ft, df01, ifelse(criteria, df1, df0))
requireNamespace("mvtnorm", quietly = TRUE)
loglik0 <- mvtnorm::dmvnorm(as.vector(object@zV - P0 %*% object@zV),
mean = rep(0, length(object@zV)), sigma = Sigma, log = TRUE)
aic0 <- -2 * loglik0 - 2 * sum(diag(P0))
bic0 <- -2 * loglik0 - sum(diag(P0)) * log(length(object@zV))
logliks <- mvtnorm::dmvnorm(as.vector(object@zV - MS %*% object@zV),
mean = rep(0, length(object@zV)), sigma = Sigma, log = TRUE)
aics <- -2 * logliks - 2 * sum(diag(SS))
bics <- -2 * logliks - sum(diag(SS)) * log(length(object@zV))
ndig <- nchar(strsplit(paste(pv),"\\.")[[1]][2])
o <- list(D = length(object@fit$psi), c = c, null = null, `tr(M0)` = sum(diag(P0)), `tr(S(phi))` = sum(diag(object@fit$S)),
`tr(S(c))` = sum(diag(SS)), `tr(M(phi))` = sum(diag(object@fit$M)), `tr(M(c))` = sum(diag(MS)),
F = round(Ft, 5), p = ifelse(ndig>6,signif(pv),round(pv, 6)), loglik0 = loglik0, aic0 = aic0, bic0 = bic0, loglik = logliks,
aic = aics, bic = bics, df0 = df0, df1 = round(df1, 5), df01 = round(df01, 5), SS0 = round(SS0, 5),
SS1 = round(SS1, 5), SS01 = round(SS01, 5))
attr(o, "SX") <- SX
attr(o, "SZ") <- SZ
attr(o, "PX") <- px
attr(o, "PZ") <- pz
attr(o, "PS") <- PS
attr(o, "A") <- Qsc
attr(o, "Sigma") <- Sigma
attr(o, "sigma") <- sqrt(object@fit$sigma2s)
attr(o, "M0") <- P0
attr(o, "P1") <- SS
attr(o, "M1") <- MS
attr(o, "B0") <- B0
attr(o, "B1") <- B1
attr(o, "B01") <- B01
attr(o, "H0") <- ifelse(null == 0, "Linear", "Spline")
attr(o, "H1") <- ifelse(null == 0, "Spline", "Linear")
return(o)
}
.scpApproxGet <- function(object, tol){
if(missing(tol)){ tol <- .Machine$double.neg.eps*1.0e-10 } else { if(is.null(tol)){ tol <- .Machine$double.neg.eps*1.0e-10 } }
cVal <- abs(sum(diag(object@fit$S))-1:nrow(object@fit$M))
cVal <- match(min(cVal),cVal)
app <- .scpLinearApprox(object = object, c = cVal, tol = tol)
att <- attributes(app)
object@fit$fitted.values <- att$M1%*%object@zV
object@fit$M <- att$M1
if(!is.null(att$PS)){
object@fit$g <- att$P1%*%(diag(1,nrow(att$P1))-att$PS)%*%object@zV
} else {
object@fit$g <- att$P1%*%object@zV
}
object@fit$approximate <- TRUE
return(object)
}
.scpTestGet <- function(object, tol){
if(missing(tol)){ tol <- .Machine$double.neg.eps*1.0e-10 } else { if(is.null(tol)){ tol <- .Machine$double.neg.eps*1.0e-10 } }
if(!is.null(object@fit$approximate))
if(object@fit$approximate) stop("Object must be a fit from scp command.")
if(is.null(object@fit$g))
stop("Model must include a spline component.")
cVal <- abs(sum(diag(object@fit$S))-1:nrow(object@fit$M))
cVal <- match(min(cVal),cVal)
app <- .scpLinearApprox(object = object, c = cVal, tol = tol)
att <- attributes(app)
app$H0 <- att$H0
app$H1 <- att$H1
app$G <- app$D
app <- app[match(c("G","c","H0","H1","df0","df1","df01","SS0","SS1","SS01","F","p"),names(app))]
tab <- cbind(c(app$df01,app$df1,app$df0),c(app$SS01,app$SS1,app$SS0),c(app$SS01/app$df01,app$SS1/app$df1,app$SS0/app$df0),c(app$F,NA,NA),c(app$p,NA,NA))
rownames(tab) <- c(ifelse(app$H1=="Spline","Approx. over Linear","Linear over Approx."),ifelse(app$H1=="Spline","Approx.","Linear"),ifelse(app$H1=="Spline","Linear","Approx."))
colnames(tab) <- c("DF", "SS", "MS", "F(y)", "Prob(F>F(y))")
print(tab,na.print="")
}
.sssFit.valid <- function(x){
el <- list(data="sss", zV="numeric", XL="ANY", XC="ANY", XF="ANY", XS="ANY", fit="list", call="call")
requireNamespace("methods",quietly=TRUE)
exists.slot <- match(names(el),methods::slotNames(x))
if(!any(is.na(exists.slot))){
if(!is.null(x@data)){
if(!is.sss(x@data)){
paste("data is not of class \"sss\".")
}
} else {
paste("a data of class \"sss\" must be part of the object.")
}
if(length(as.vector(unlist(x@zV)))!=nrow(x@data@coords)){
paste("Wrong structure in zV and data.")
}
if(TRUE){
valid = NULL
for(i in 1:length(el)){
if(!is.na(match(el[[i]],c("sss","numeric","list","call")))){
valid[i] <- any(grepl(paste("^",el[[i]],"$",sep=""),class(methods::slot(x,names(el)[[i]]))))
} else {
valid[i] = TRUE
}
}
if(all(valid)) TRUE else paste("Slots",paste(names(el)[!valid],sep="",collapse=", "),"not match their class.")
}
} else {
paste("Slots",paste(names(el)[is.na(exists.slot)],sep="",collapse=", "),"are missing.")
}
}
sssFit <- setClass("sssFit",
slots = list(data="sss", zV="numeric", XL="ANY", XC="ANY", XF="ANY", XS="ANY", fit="list", call="call"),
package = "scpm", validity = .sssFit.valid)
.sssToSurface = function (data = NULL, grid.list = NULL, values = NULL, names = NULL){
if(!is.null(data)){
grid.list = attr(data@grid, "grid.list")
}
if(!is.null(grid.list)){
if(is.matrix(grid.list)){
grid.list = attr(grid.list, "grid.list")
}
}
if(is.null(names)){
if(!is.null(values)){
list(x = grid.list[[1]], y = grid.list[[2]], z = matrix(values, ncol = length(grid.list[[2]]), nrow = length(grid.list[[1]]), byrow = FALSE), xlab = names(grid.list)[1], ylab = names(grid.list)[2])
} else {stop("values or names must be specified!")}
} else {
out = NULL
for(i in 1:length(names)){
out[[i]] = list(x = grid.list[[1]], y = grid.list[[2]], z = matrix(as.vector(data@data[[names[i]]]), ncol = length(grid.list[[2]]), nrow = length(grid.list[[1]]), byrow = FALSE), xlab = names(grid.list)[1], ylab = names(grid.list)[2])
}
names(out) = names
return(out)
}
}
.add3Dpoints <- function(x,y,z, pmat) {
tr <- cbind(x,y,z,1) %*% pmat
list(x = tr[,1]/tr[,4], y= tr[,2]/tr[,4])
}
.legendMapSimple <- function(col, lev, val=NULL){
requireNamespace("graphics", quietly = TRUE)
opar <- graphics::par
n <- length(col)
range.lev = c(min(lev),max(lev))
bx <- graphics::par("usr")
box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000,
bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
box.cy <- c(bx[3], bx[3])
box.sy <- (bx[4] - bx[3]) / n
unit.lev = (bx[4] - bx[3]) / (range.lev[2] - range.lev[1])
olev = lev[order(lev)]
xx <- rep(box.cx, each = 2)
xx = xx - (xx[4]-xx[1])/2
graphics::par(xpd = TRUE)
yb = NULL
for(i in 1:n){
yb[i] = box.cy[1] + box.sy/2
yy <- c(box.cy[1] + (box.sy * (i - 1)),
box.cy[1] + (box.sy * (i)),
box.cy[1] + (box.sy * (i)),
box.cy[1] + (box.sy * (i - 1)))
yy = yy + box.sy/2
yb[i] = mean(yy)
graphics::polygon(xx, yy, col = col[i], border = col[i])
}
requireNamespace("stats", quietly = TRUE)
lev.tick = stats::quantile(lev, probs = seq(0,1,0.1))
lev.plot = lev[match(paste(lev.tick),paste(lev))]
yb.plot = yb[match(paste(lev.tick),paste(lev))]
valb.plot = min(yb) + (val-min(val))*(max(yb)-min(yb))/(max(val)-min(val))
graphics::axis(side = 4, at = yb.plot, labels = round(lev.plot), las = 3, tick = TRUE, line = 1, tcl = -0.5, cex.axis = 0.85)
graphics::axis(side = 4, at = valb.plot, labels = FALSE, las = NULL, tick = TRUE, line = 1, lwd.ticks = 2, tcl = -0.25)
par <- opar
}
.legendMap <- function(col, lev, val=NULL, by.prob = 0.01, out = FALSE){
requireNamespace("graphics", quietly = TRUE)
opar <- graphics::par
n <- length(col)
range.lev = c(min(lev, na.rm=TRUE),max(lev, na.rm=TRUE))
bx <- graphics::par("usr")
box.cx <- c(bx[2] + (bx[2] - bx[1]) / 1000,
bx[2] + (bx[2] - bx[1]) / 1000 + (bx[2] - bx[1]) / 50)
box.cy <- c(bx[3], bx[3])
box.sy <- (bx[4] - bx[3]) / n
unit.lev = (bx[4] - bx[3]) / (range.lev[2] - range.lev[1])
olev = lev[order(lev)]
xx <- rep(box.cx, each = 2)
xx = xx - (xx[4]-xx[1])/2
graphics::par(xpd = TRUE)
yb = NULL
for(i in 1:n){
yb[i] = box.cy[1] + box.sy/2
yy <- c(box.cy[1] + (box.sy * (i - 1)),
box.cy[1] + (box.sy * (i)),
box.cy[1] + (box.sy * (i)),
box.cy[1] + (box.sy * (i - 1)))
yy = yy + box.sy/2
yb[i] = mean(yy)
graphics::polygon(xx, yy, col = col[i], border = col[i])
}
requireNamespace("stats", quietly = TRUE)
lev.tick = stats::quantile(lev, probs = seq(0,1,by.prob), na.rm = TRUE)
lev.plot = lev[match(paste(lev.tick),paste(lev))]
yb.plot = yb[match(paste(lev.tick),paste(lev))]
if(!is.null(val)){
valb.plot = min(yb) + (val-min(val))*(max(yb)-min(yb))/(max(val)-min(val))
}
graphics::axis(side = 4, at = yb.plot, labels = round(lev.plot), las = 3, tick = TRUE, line = 1, tcl = -0.5, cex.axis = 0.85)
if(!is.null(val)){
graphics::axis(side = 4, at = valb.plot, labels = FALSE, las = NULL, tick = TRUE, line = 1, lwd.ticks = 2, tcl = -0.25)
}
par <- opar
if(out){
list(usr = yb)
}
}
.colorVecMat <- function(obj, which = c("rainbow","heat.colors","terrain.colors","topo.colors","cm.colors","colorRamp","colorRampPalette"), col.args = list()){
.defaultArg("which","rainbow")
which.fx = c("colorRamp","colorRampPalette")
if(is.matrix(obj)){
values = as.vector(obj)
} else {
values = obj
}
requireNamespace("grDevices", quietly = TRUE)
cargs = list()
if(length(col.args)>0){
for(i in 1:length(col.args)){
cargs[[i]] = col.args[[i]]
}
names(cargs) = names(col.args)
}
if(!is.na(match(which,which.fx))){
af = formals(which)
fx.args = cargs[!is.na(match(names(cargs),names(af)))]
fx.cols = do.call(which,fx.args)
which = "fx.cols"
}
if(is.null(cargs$n)){
cargs$n = length(values)
}
af = formals(which)
cargs = cargs[!is.na(match(names(cargs),names(af)))]
cols = do.call(which,cargs)
if(is.matrix(obj)){
zfacet <- obj[-1, -1] + obj[-1, -ncol(obj)] + obj[-nrow(obj), -1] + obj[-nrow(obj), -ncol(obj)]
facetcol <- cut(zfacet, length(cols))
return(cols[facetcol])
} else {
oz = values[order(values)]
icols = match(paste(values),paste(oz))
return(cols[icols])
}
}
.scpImage <- function(values, data = NULL, grid.list = NULL,
which = c("image","image.default","image.plot","image.kriging"),
which.col = "rainbow", image.args = list(), contour.args = list(),
legend.args = list(), col.args = list()){
.defaultArg("which","image.default")
if(is.list(values)){
zsurf = values
} else {
zsurf = .sssToSurface(data = data, grid.list = grid.list, values = values, names = NULL)
}
v.v = as.vector(zsurf$z)
v.m = zsurf$z
if(is.null(col.args$start)){ col.args$start = min(.transfTo01(v.v)) }
if(is.null(col.args$end)){ col.args$end = max(.transfTo01(v.v)) }
if(is.null(image.args$col)){ image.args$col = .colorVecMat(obj = v.v[order(v.v)], which = which.col, col.args = col.args) }
if(is.null(image.args$oldstyle)){ image.args$oldstyle = TRUE }
if(is.null(image.args$useRaster)){ image.args$useRaster = TRUE }
if(is.null(contour.args$method)){ contour.args$method = "flattest" }
if(is.null(legend.args$col)){ legend.args$col = image.args$col }
if(is.null(legend.args$lev)){ legend.args$lev = v.v[order(v.v)] }
if(is.null(legend.args$val)){ legend.args$val = v.v[order(v.v)] }
if(is.null(legend.args$by.prob)){ legend.args$by.prob = 0.01 }
if(is.null(legend.args$out)){ legend.args$out = FALSE }
af = formals(which)
image.args = image.args[!is.na(match(names(image.args),names(af)))]
af = formals(contour.default)
contour.args = contour.args[!is.na(match(names(contour.args),names(af)))]
iargs = zsurf
cargs = zsurf
for(i in 1:length(image.args)){
iargs[[i+length(zsurf)]] = image.args[[i]]
}
names(iargs) = c(names(zsurf),names(image.args))
if(grepl("^image.plot$",which))
if(!requireNamespace("fields", quietly = TRUE)) which <- "image.default"
if(grepl("^image.kriging$",which))
which <- "image.default"
if(grepl("^image(.default){0,1}$",which))
requireNamespace("graphics", quietly = TRUE)
do.call(which,iargs)
for(i in 1:length(contour.args)){
cargs[[i+length(zsurf)]] = contour.args[[i]]
}
names(cargs) = c(names(zsurf),names(contour.args))
cargs$add = TRUE
if(!isNamespaceLoaded("graphics"))
requireNamespace("graphics", quietly = TRUE)
do.call("contour.default",cargs)
if(which!="image.plot"){
do.call(".legendMap",legend.args)
}
}
.scpPersp = function(values, data = NULL, grid.list = NULL,
which = c("persp","persp.kriging","persp.grf"),
which.col = "rainbow", persp.args = list(), legend.args = list(), col.args = list()){
.defaultArg("which","persp")
if(is.list(values)){
zsurf = values
} else {
zsurf = .sssToSurface(data = data, grid.list = grid.list, values = values, names = NULL)
}
v.v = as.vector(zsurf$z)
v.m = zsurf$z
if(is.null(col.args$start)){ col.args$start = min(.transfTo01(v.v)) }
if(is.null(col.args$end)){ col.args$end = max(.transfTo01(v.v)) }
if(is.null(persp.args$col)){ persp.args$col = .colorVecMat(v.m, which = which.col, col.args = col.args) }
if(is.null(persp.args$theta)){ persp.args$theta = 0 }
if(is.null(persp.args$phi)){ persp.args$phi = 60 }
if(is.null(persp.args$ticktype)){ persp.args$ticktype = "detailed" }
if(is.null(persp.args$border)){ persp.args$border = "transparent" }
if(is.null(persp.args$shade)){ persp.args$shade = 0.15 }
if(is.null(persp.args$zlab)){ persp.args$zlab = "z" }
if(is.null(legend.args$col)){ legend.args$col = .colorVecMat(v.v[order(v.v)], which = which.col,
col.args = col.args) }
if(is.null(legend.args$lev)){ legend.args$lev = v.v[order(v.v)] }
if(is.null(legend.args$val)){ legend.args$val = v.v[order(v.v)] }
if(is.null(legend.args$by.prob)){ legend.args$by.prob = 0.01 }
if(is.null(legend.args$out)){ legend.args$out = FALSE }
pargs = zsurf
for(i in 1:length(persp.args)){
pargs[[i+length(zsurf)]] = persp.args[[i]]
}
names(pargs) = c(names(zsurf),names(persp.args))
if(grepl("^persp.(grf|kriging){1}$",which))
which <- "persp"
if(grepl("^persp$",which))
requireNamespace("graphics", quietly = TRUE)
xypos = do.call(which,pargs)
do.call(".legendMap",legend.args)
return(xypos)
}
.scpPoints = function(obj, values = NULL, pos3d = NULL, points.args = list()){
if(is.sss(obj)){
coords = obj@coords
}
if(is.matrix(obj)){
coords = obj
}
af = formals(points)
points.args = points.args[!is.na(match(names(points.args),names(af)))]
pargs = list()
if(is.null(values) & is.null(pos3d)){
pargs$x = coords[,1]
pargs$y = coords[,2]
} else {
pargs = .add3Dpoints(coords[,1],coords[,2], values, pos3d)
}
if(length(points.args)>0){
bar = length(pargs)
nbar = names(pargs)
for(i in 1:length(points.args)){
pargs[[i+bar]] = points.args[[i]]
}
names(pargs) = c(nbar,names(points.args))
}
requireNamespace("graphics", quietly = TRUE)
do.call("graphics::points",pargs)
}
.scpPlot = function(obj, data = NULL,
which.col = "rainbow", which.image = "image.default", which.persp = "persp",
image.args = list(), contour.args = list(), persp.args = list(), col.args = list(),
add.points = FALSE,
titles = list()){
if(is.null(data)){
data = obj$data
}
if(is.null(titles$z)){ titles$z = bquote(Observed~Field) }
if(is.null(titles$mu)){ titles$mu = bquote(Estimated~Field,~E(Z)==mu[Z]) }
if(is.null(titles$g)){ titles$g = bquote(Estimated~g(x,y)) }
if(is.null(titles$cb)){ titles$cb = bquote(Estimated~Cb) }
if(is.null(titles$xb)){ titles$xb = bquote(Estimated~X*beta) }
add.one = 0
c.mu = obj$fit$fitted.values
if(ncol(obj$fit$X)>4){
c.cb = obj$fit$X[,1:(ncol(obj$fit$X)-4)]%*%obj$fit$beta[1:(ncol(obj$fit$X)-4)]
add.one = add.one + 1
} else {
c.cb = NULL
}
c.gf = obj$fit$g
c.xb = obj$fit$X[,(ncol(obj$fit$X)-3):ncol(obj$fit$X)]%*%obj$fit$beta[(ncol(obj$fit$X)-3):ncol(obj$fit$X)]
c.all = c(obj$zV, c.mu, c.cb, c.gf, c.xb)
out.z = .sssToSurface(data = data, values = obj$zV)
out.mu = .sssToSurface(data = data, values = c.mu)
if(ncol(obj$fit$X)>4){
out.cb = .sssToSurface(data = data, values = c.cb)
}
out.gf = .sssToSurface(data = data, values = c.gf)
out.xb = .sssToSurface(data = data, values = c.xb)
par(mfrow=c(2,ifelse(add.one>0,5,3)), mar=c(5,4,4,3)+0.1)
.scpImage(values = out.z, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$z)
if(add.one>0){
.scpImage(values = out.mu, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$mu)
.scpImage(values = out.gf, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$g)
.scpImage(values = out.cb, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$cb)
.scpImage(values = out.xb, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$xb)
} else {
.scpImage(values = out.gf, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$g)
.scpImage(values = out.xb, which = which.image, which.col = which.col,
image.args = image.args, contour.args = contour.args, col.args = col.args)
title(main = titles$xb)
}
if(is.null(persp.args$zlim)){ persp.args$zlim = c(min(c.all),max(c.all)) }
pp = .scpPersp(values = out.z, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$z)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
if(add.one>0){
pp = .scpPersp(values = out.mu, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$mu)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
pp = .scpPersp(values = out.gf, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$g)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
pp = .scpPersp(values = out.cb, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$cb)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
pp = .scpPersp(values = out.xb, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$xb)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
} else {
pp = .scpPersp(values = out.gf, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$g)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
pp = .scpPersp(values = out.xb, which = which.persp, which.col = which.col, persp.args = persp.args, col.args = col.args)
title(main = titles$xb)
if(add.points){ .scpPoints(obj = data, values = obj$zV, pos3d = pp, points.args = list(pch = 1)) }
}
}
.getSemivariogram <- function(x,obj){
.semiVarF(h=x,phi=obj@fit$phi,kappa=obj@fit$kappa,model=eval(obj@call$model),sill=obj@fit$sigma2,nugget=obj@fit$tau2,tol.nugget=ifelse(is.null(eval(obj@call$nugget.tol)),1.0e-15,eval(obj@call$nugget.tol)),use.cor=TRUE)
}
.plotSemivariogram <- function(obj,h=NULL,hmax=NULL,size=200,plot=TRUE,prange=TRUE,legend=TRUE,...){
aux.tol <- ifelse(is.null(eval(obj@call$nugget.tol)),1.0e-15,eval(obj@call$nugget.tol))
if(is.null(hmax) & is.null(h))
hmax <- max(abs(obj@XS$h)) + (max(abs(obj@XS$h))-min(abs(obj@XS$h)))*0.05
if(is.null(hmax) & !is.null(h))
hmax <- max(h)
if(is.null(h)){
xv <- seq(0,hmax,l=size)
} else {
xv <- h
}
sv <- .getSemivariogram(xv,obj)
xv <- xv[!is.na(sv)]
sv <- sv[!is.na(sv)]
xv <- c(aux.tol,xv)
sv <- c(obj@fit$tau2,sv)
sv <- sv[order(xv)]
xv <- xv[order(xv)]
rr <- xv[match(max(sv),sv)]
rr <- (xv[sv==as.numeric(obj@fit$sigma2s)])[1]
rp <- xv[sv <= as.numeric(obj@fit$sigma2s)*0.95]
rp <- rp[length(rp)]
if(plot){
message(paste("Sill=",round(obj@fit$sigma2s,6),sep=""))
message(paste("Practical range h",ifelse(length(rp)==0,paste(">",hmax,sep=""),paste("=",round(rp,6),sep="")),", 0.95*sill=",round(as.numeric(obj@fit$sigma2s)*0.95,6),sep=""))
gargs <- list(...)
if(is.null(gargs$xlim)){xlim <- c(0,max(xv))} else {xlim <- gargs$xlim}
if(is.null(gargs$ylim)){ylim <- c(0,max(sv))} else {ylim <- gargs$ylim}
if(is.null(gargs$main)){main <- paste("Fitted semi-variogram model",eval(obj@call$model))} else {main <- gargs$main}
plot(sv[xv>=aux.tol]~xv[xv>=aux.tol],ylab=expression(hat(gamma(h))),xlab=expression(h),type="l",xlim=xlim,ylim=ylim,main=main)
points(x=aux.tol,y=obj@fit$tau2,pch=16)
if(prange){
rr <- NA
abline(v=c(if(!is.na(rr)){rr},if(length(rp)!=0){rp}),lty=c(if(!is.na(rr)){4},if(length(rp)!=0){3}),col=c(if(!is.na(rr)){"blue"},if(length(rp)!=0){"blue"}))
if(legend)
legend("bottomright",legend=c(if(!is.na(rr)){as.expression(bquote(Range==.(round(rr,4))))},if(length(rp)!=0){as.expression(bquote(Prac.*phantom(0)*Range==.(round(rp,4))))}),lty=c(if(!is.na(rr)){4},if(length(rp)!=0){3}),col=c(if(!is.na("blue")){rr},if(length(rp)!=0){"blue"}),bg="white")
}
} else {
list(x = xv[xv>=aux.tol], y = sv[xv>=aux.tol], sill = obj@fit$sigma2s, prange = rp, psill = obj@fit$sigma2s*0.95)
}
}
.scpSummary <- function(object,alpha=0.05){
if(!is.null(object@fit$approximate))
if(object@fit$approximate) stop("Object must be a fit from scp command.")
obj <- object
tv <- obj@fit$beta/sqrt(diag(obj@fit$var.beta))
ptv <- 2*(1-pt(abs(tv),length(obj@zV)-ncol(obj@fit$X)))
atv <- qt(1-alpha/2,length(obj@zV)-ncol(obj@fit$X))*sqrt(diag(obj@fit$var.beta))
ctable <- cbind(obj@fit$beta,sqrt(diag(obj@fit$var.beta)),tv,ptv,obj@fit$beta-atv,obj@fit$beta+atv)
colnames(ctable) <- c("estimate","std.error","t(y)","p(|t(y)|>t)","LL(95%)","UL(95%)")
nr <- rownames(ctable)
nr <- gsub("^Intercept$","(Intercept)",nr)
rownames(ctable) <- nr
return(ctable)
}
.scpImagePlot <- function(obj, type = c("obs","fit","g","both"),...){
if(missing(type)){ type <- "fit"} else { if(is.null(type)){ type <- "fit" } }
type <- match.arg(type)
requireNamespace("interp", quietly = TRUE)
intz <- interp::interp(x=obj@data@coords[,1], y=obj@data@coords[,2], z=obj@zV, method = "linear", linear = TRUE, extrap = FALSE)
intm <- interp::interp(x=obj@data@coords[,1], y=obj@data@coords[,2], z=obj@fit$fitted.values, method = "linear", linear = TRUE, extrap = FALSE)
if(!is.null(obj@fit$g))
intg <- interp::interp(x=obj@data@coords[,1], y=obj@data@coords[,2], z=obj@fit$g, method = "linear", linear = TRUE, extrap = FALSE)
gargs <- list(...)
if(type=="both")
par(mfrow=c(1,2))
if(type=="obs" | type=="both"){
image(intz,...)
if(is.null(gargs$main)){
title("Observed values")
}
contour(intz,add=TRUE)
}
if(type=="g" | type=="both"){
if(!is.null(obj@fit$g)){
image(intg,...)
if(is.null(gargs$main)){
title("Fitted spline")
}
contour(intg,add=TRUE)
}
}
if(type=="fit" | type=="both"){
image(intm,...)
if(is.null(gargs$main)){
title("Fitted model")
}
contour(intm,add=TRUE)
}
}
.scpLevelPlot <- function(obj,what=c("obs","fit","g"),level.at="fivenum",colors=c("yellow","red"),...){
if(missing(what)){ what <- "fit"} else { if(is.null(what)){ what <- "fit" } }
what <- match.arg(what)
requireNamespace("interp",quietly=TRUE)
if(what=="obs")
values <- obj@zV
if(what=="fit")
values <- obj@fit$fitted.values
if(what=="g")
values <- obj@fit$g
int <- interp::interp(x=obj@data@coords[,1], y=obj@data@coords[,2], z=values, method = "linear", linear = TRUE, extrap = FALSE)
igridxy <- (cbind(1,int$y)%x%cbind(1,int$x))[,2:3]
igridx <- igridxy[,1]
igridy <- igridxy[,2]
igridz <- as.vector(int$z)
if(!is.null(level.at)){
if(is.character(level.at)){
level.at <- as.vector(apply(cbind(igridz),2,level.at))
}
} else {
level.at <- fivenum(igridz)
}
lev.color <- .colorVecMat(level.at, which = "colorRampPalette", col.args = list(colors=colors))
requireNamespace("lattice",quietly=TRUE)
lattice::levelplot(igridz~igridx*igridy, at = level.at, col.regions = lev.color, panel = lattice::panel.levelplot.raster, xlab = names(obj@data@coords)[1], ylab = names(obj@data@coords)[2], ...)
}
.scpModelPlot <- function(x, what = c("obs","fit","g"), type = c("levelplot", "image", "persp", "persp3d"), which = "colorRampPalette", col.args = list(colors = c("yellow","red")), col.contour = "black", level.at = "fivenum", border = "transparent", theta = 0, phi = 45, shade = 0.1, ...){
obj <- x
if(missing(what)){ what <- "fit" } else { if(is.null(what)) { what <- "fit" }}
what <- match.arg(what)
if(missing(type)){ type <- "image" } else { if(is.null(type)) { type <- "image" }}
type <- match.arg(type)
if(type=="levelplot"){
.scpLevelPlot(obj, level.at = level.at, what = what)
} else {
requireNamespace("interp", quietly = TRUE)
requireNamespace("lattice", quietly = TRUE)
if(what=="obs")
vals <- obj@zV
if(what=="fit")
vals <- obj@fit$fitted.values
if(what=="g"){
if(!is.null(obj@fit$g)){
vals <- obj@fit$g
} else {
return(message("No g function in the model!"))
}
}
int <- interp::interp(x=obj@data@coords[,1], y=obj@data@coords[,2], z=vals, method = "linear", linear = TRUE, extrap = FALSE)
vv <- int$z
if(type=="image" | type=="levelplot"){
vv <- as.vector(vv)[order(vv)]
cols <- .colorVecMat(vv, which = which, col.args = col.args)
} else if(type=="persp"){
cols <- .colorVecMat(vv, which = which, col.args = col.args)
} else if(type=="persp3d"){
vv <- as.vector(vv)
cols <- .colorVecMat(vv, which = which, col.args = col.args)
}
cf <- match.call(expand.dots=TRUE)
cto <- cf[-c(na.omit(match(c("obj","what","type","which","col.args","col.contour","level.at",if(type!="persp" & type!="persp3d"){c("border")},if(type!="persp"){c("theta","phi","shade")}),names(cf))))]
cto$x <- int$x
cto$y <- int$y
cto$z <- int$z
cto$col <- cols
if(is.null(cto$xlab))
cto$xlab <- names(obj@data@coords)[1]
if(is.null(cto$ylab))
cto$ylab <- names(obj@data@coords)[2]
if(is.null(cto$zlab) & type!="image" & type!="levelplot")
cto$zlab <- "z"
if(type=="image"){
cto[[1]] <- quote(image)
eval(cto)
cto[[1]] <- quote(contour)
cto$col <- col.contour
if(is.null(cto$add))
cto$add <- TRUE
eval(cto)
}
if(type=="persp"){
if(all(na.omit(cto$z==mean(cto$z,na.rm=TRUE))))
cto$zlim <- c(min(cto$z,na.rm=TRUE)-sqrt(.Machine$double.neg.eps),max(cto$z,na.rm=TRUE)+sqrt(.Machine$double.neg.eps))
cto[[1]] <- quote(persp)
eval(cto)
}
if(type=="persp3d"){
requireNamespace("rgl",quietly=TRUE)
cto[[1]] <- quote(rgl::persp3d)
eval(cto)
}
}
}
.getInformationCriterion <- function(object,what=c("m-aic","c-aic","j-bic","gic","hq-gic","b-gic","pn-gic","aic","bic","c-bic"),k=2,tol=.Machine$double.neg.eps){
obj <- object
if(!is.null(object@fit$approximate))
if(object@fit$approximate) stop("Object must be a fit from scp command.")
if(missing(what)){ what <- "aic"} else { if(is.null(what)){ what <- "aic" } }
what <- match.arg(what)
np <- ncol(obj@fit$X)
nep <- sum(diag(obj@fit$M))
nq <- length(c(obj@fit$sigma2s,if(!is.null(obj@fit$alpha)){obj@fit$alpha},if(!is.numeric(eval(obj@call$fix.nugget))){obj@fit$tau2},obj@fit$phi,if(!is.numeric(eval(obj@call$fix.kappa))){obj@fit$kappa}))
Sigma <- obj@fit$sigma2s%x%obj@fit$R
mu <- obj@fit$fitted.values
Xb <- mu
V <- Sigma
if(!is.null(obj@fit$Zs) & !is.null(obj@fit$r)){
Xb <- mu - obj@fit$Zs%*%obj@fit$r
V <- Sigma + obj@fit$sigma2s%x%obj@fit$Zs%*%obj@fit$penalty$decompose$A%*%t(obj@fit$Zs)
}
if(what=="m-aic"){#Vaida and Blanchard (2005)
a0 <- k*length(obj@zV)/(length(obj@zV)-np-nq-1)
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="c-aic"){#Burnham and Anderson (2002)
a0 <- k
npar <- nep + nq
vmean <- mu
mvar <- Sigma
}
if(what=="aic"){#Akaike's
a0 <- k
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="bic"){#Schwart's
a0 <- log(length(obj@zV))
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="c-bic"){#Mario's
a0 <- log(length(obj@zV))
npar <- nep + nq
vmean <- mu
mvar <- Sigma
}
if(what=="j-bic"){#Jones (2011)
vone <- rep(1,length(obj@zV))
U <- diag(diag(V),nrow(V))
a0 <- log(t(vone)%*%U%*%solve(V,tol=tol)%*%U%*%vone)
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="gic"){#Pu and Niu (2006)
a0 <- k
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="hq-gic"){#Hannan-Quinn (1979)
a0 <- k*log(log(length(obj@zV)))
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="b-gic"){#Bozgodan (1987)
a0 <- (log(length(obj@zV))+1)
npar <- np + nq
vmean <- Xb
mvar <- V
}
if(what=="pn-gic"){#Pu and Niu (2006)
a0 <- sqrt(log(length(obj@zV)))
npar <- np + nq
vmean <- Xb
mvar <- V
}
pen <- a0*npar
requireNamespace("mvtnorm",quietly=TRUE)
lL <- mvtnorm::dmvnorm(x=obj@zV,mean=vmean,sigma=mvar,log=TRUE)
cr <- 2*lL - pen
list(logLik = as.numeric(lL), criterion = as.numeric(-cr), ka0 = as.numeric(a0), numpar = npar, penalty = as.numeric(pen))
}
setGeneric("scpApproximate",function(object, tol){ standardGeneric("scpApproximate") })
setMethod("scpApproximate", signature(object = "sssFit"), function(object, tol){
.scpApproxGet(object = object, tol = tol)
})
setGeneric("testSurface",function(object, tol){ standardGeneric("testSurface") })
setMethod("testSurface", signature(object = "sssFit"), function(object, tol){
.scpTestGet(object = object, tol = tol)
})
setGeneric("plot")
setGeneric("summary")
setGeneric("AIC")
setGeneric("BIC")
setGeneric("AICm",function(object, k, only.criterion){ standardGeneric("AICm") })
setGeneric("AICc",function(object, k, only.criterion){ standardGeneric("AICc") })
setGeneric("BICc",function(object, only.criterion){ standardGeneric("BICc") })
setGeneric("BICj",function(object, k, tol, only.criterion){ standardGeneric("BICj") })
setGeneric("GIC",function(object, k, only.criterion){ standardGeneric("GIC") })
setGeneric("GIChq",function(object, k, only.criterion){ standardGeneric("GIChq") })
setGeneric("GICpn",function(object, only.criterion){ standardGeneric("GICpn") })
setGeneric("GICb",function(object, only.criterion){ standardGeneric("GICb") })
setGeneric("Variogram",function(object, distance, plot, ...){ standardGeneric("Variogram") })
setMethod("plot", signature(x = "sssFit", y = "missing"), function(x, what, type, which, col.args, col.contour, level.at, border, theta, phi, shade, ...){
if(missing(what)){ what <- "fit" } else {if(is.null(what)){ what <- "fit" }}
if(missing(type)){ type <- "image" } else {if(is.null(type)){ type <- "image" }}
if(missing(which)){ which <- "colorRampPalette" } else {if(is.null(which)){ which <- "colorRampPalette" }}
if(missing(col.args)){ col.args <- list(colors=c("yellow","red")) } else {if(is.null(col.args)){ col.args <- list(colors=c("yellow","red")) }}
if(missing(col.contour)){ col.contour <- "black" } else {if(is.null(col.contour)){ col.contour <- "black" }}
if(missing(level.at)){ level.at <- "fivenum" } else {if(is.null(level.at)){ level.at <- "fivenum" }}
if(missing(border)){ border <- "transparent" } else {if(is.null(border)){ border <- "transparent" }}
if(missing(theta)){ theta <- 0 } else {if(is.null(theta)){ theta <- 0 }}
if(missing(phi)){ phi <- 45 } else {if(is.null(phi)){ phi <- 45 }}
if(missing(shade)){ shade <- 0.1 } else {if(is.null(shade)){ shade <- 0.1 }}
.scpModelPlot(x, what, type, which, col.args, col.contour, level.at, border, theta, phi, shade, ...)
})
setMethod("summary", signature(object = "sssFit"), .scpSummary)
setMethod("AIC", signature(object = "sssFit", k = "ANY"), function(object, k, only.criterion){
if(missing(k)){ k = 2 } else { if(is.null(k)){ k = 2 } }
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="aic",k=k)$criterion
} else {
.getInformationCriterion(object,what="aic",k=k)
}
})
setMethod("BIC", signature(object = "sssFit"), function(object, only.criterion){
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="bic")$criterion
} else {
.getInformationCriterion(object,what="bic")
}
})
setMethod("AICm", signature(object = "sssFit", k = "ANY"), function(object, k, only.criterion){
if(missing(k)){ k = 2 } else { if(is.null(k)){ k = 2 } }
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="m-aic",k=k)$criterion
} else {
.getInformationCriterion(object,what="m-aic",k=k)
}
})
setMethod("AICc", signature(object = "sssFit", k = "ANY"), function(object, k, only.criterion){
if(missing(k)){ k = 2 } else { if(is.null(k)){ k = 2 } }
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="c-aic",k=k)$criterion
} else {
.getInformationCriterion(object,what="c-aic",k=k)
}
})
setMethod("BICc", signature(object = "sssFit"), function(object, only.criterion){
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="c-bic")$criterion
} else {
.getInformationCriterion(object,what="c-bic")
}
})
setMethod("BICj", signature(object = "sssFit", k = "ANY", tol = "ANY"), function(object, k, tol, only.criterion){
if(missing(k)){ k = 2 } else { if(is.null(k)){ k = 2 } }
if(missing(tol)){ tol = .Machine$double.neg.eps } else { if(is.null(tol)){ tol = .Machine$double.neg.eps } }
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="j-bic",k=k,tol=tol)$criterion
} else {
.getInformationCriterion(object,what="j-bic",k=k,tol=tol)
}
})
setMethod("GIC", signature(object = "sssFit", k = "ANY"), function(object, k, only.criterion){
if(missing(k)){ k = 2 } else { if(is.null(k)){ k = 2 } }
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="gic",k=k)$criterion
} else {
.getInformationCriterion(object,what="gic",k=k)
}
})
setMethod("GIChq", signature(object = "sssFit", k = "ANY"), function(object, k, only.criterion){
if(missing(k)){ k = 2 } else { if(is.null(k)){ k = 2 } }
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="hq-gic",k=k)$criterion
} else {
.getInformationCriterion(object,what="hq-gic",k=k)
}
})
setMethod("GICpn", signature(object = "sssFit"), function(object, only.criterion){
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="pn-gic")$criterion
} else {
.getInformationCriterion(object,what="pn-gic")
}
})
setMethod("GICb", signature(object = "sssFit"), function(object, only.criterion){
if(missing(only.criterion)){ only.criterion = TRUE } else { if(is.null(only.criterion)){ only.criterion = TRUE } }
if(only.criterion){
.getInformationCriterion(object,what="b-gic")$criterion
} else {
.getInformationCriterion(object,what="b-gic")
}
})
setMethod("Variogram", signature(object = "sssFit", distance = "ANY"),
function(object, distance, plot, ...){
if(missing(distance)){ distance <- NULL } else { if(!is.null(distance)){ if(!is.numeric(distance)){
warning("Variogram: distance must be numeric or null. Setting to NULL.")
distance <- NULL
} } }
if(missing(plot)){ plot <- TRUE } else { if(is.null(TRUE)){ plot <- TRUE } }
.plotSemivariogram(object,h=distance,hmax=NULL,size=10000,plot=plot,prange=FALSE,legend=FALSE, ...)
})
.tryCatchWE = function (expr){
W <- NULL
w.handler <- function(w) {
W <<- w
invokeRestart("muffleWarning")
}
value = withCallingHandlers(tryCatch(expr, error = function(e) e), warning = w.handler)
isE = !is.na(pmatch("Error",paste(value,collapse="")))
list(value = value, error = isE, warning = W)
}
.defaultArg = function(arg,default){
env = parent.frame()
if(!is.character(arg)) stop("arg must be character!")
txt = NULL
txt[1] = paste("if(missing(",arg,") || is.null(",arg,")){",arg," = ",shQuote(default),"}",sep="")
txt[2] = paste("if(!is.character(",arg,")){",arg," = as.character(",arg,")}",sep="")
txt[3] = paste(arg," = match.arg(",arg,")",sep="")
txt[4] = paste("if(!is.character(",default,")){",arg," = as.numeric(",arg,")}",sep="")
txt = txt[c(1:ifelse(is.character(default),3,4))]
eval(parse(text=txt),envir=env)
}
.invSchur = function(x, Re.values = FALSE, ...){
requireNamespace("Matrix", quietly = TRUE)
ex = Matrix::Schur(x, ...)
if(any(is.complex(ex$EValues))){
warning("Complex eigenvalues!")
}
if(Re.values){
d.v = Re(ex$EValues)
ex.Q = Re(ex$Q)
} else {
d.v = ex$EValues
ex.Q = ex$Q
}
ex.Q%*%diag(1/d.v)%*%t(ex.Q)
}
.invEigen = function(x, Re.values = FALSE, ...){
ex = eigen(x, ...)
if(Re.values){
d.v = Re(ex$values)
ex.v = Re(ex$vectors)
} else {
d.v = ex$values
ex.v = ex$vectors
}
ex.v%*%diag(1/d.v)%*%t(ex.v)
}
.invGsvd = function(x, tol = sqrt(.Machine$double.eps), ...){
ex = svd(x, ...)
r.s = ex$d > tol*ex$d[1]
if(any(r.s)){
ex$v[,r.s]%*%(t(ex$u[,r.s])/ex$d[r.s])
} else {x}
}
.whichInverse = function(name = c("solve","ginv",".invEigen",".invSchur",".invGsvd"), ...){
.defaultArg("name","solve")
if(name=="ginv"){
requireNamespace("MASS", quietly = TRUE)
name = paste("MASS::",name,sep="")
}
l.args = list(...)
o.args = formals(get("name"))
c.args = l.args[-1]
v.args = match(names(c.args),names(o.args))
if(any(!is.na(v.args))){
l.args = l.args[c(TRUE,!is.na(v.args))]
} else {
l.args = l.args[1]
}
do.call(name, l.args)
}
.transfTo01 <- function(y,dom=NULL,to.Inf=.Machine$double.base**(-.Machine$double.min.exp)){
dom.obs <- c(min(y),max(y))
if(is.null(dom)){
dom <- dom.obs
}
lim.a <- ifelse(min(dom) <= -to.Inf, -to.Inf, min(dom))
lim.b <- ifelse(max(dom) >= +to.Inf, +to.Inf, max(dom))
y.01 <- (y - lim.a)/(lim.b-lim.a)
return(y.01)
}
.transfFrom01 <- function(y,dom,to.Inf=.Machine$double.base**(-.Machine$double.min.exp)){
if(missing(dom))
stop("domain of the original values must be specified (otherwise infinity transformations are available)")
lim.a <- ifelse(min(dom) <= -to.Inf, -to.Inf, min(dom))
lim.b <- ifelse(max(dom) >= +to.Inf, +to.Inf, max(dom))
y.or <- y*(lim.b-lim.a) + lim.a
return(y.or)
}
.getKnots <- function(x){
requireNamespace("stats", quietly = TRUE)
keep.class = class(x)
aux = stats::xtabs(~x)
knots.values <- names(aux)
class(knots.values) <- keep.class
return(knots.values)
}
.tpsBase = function(x,...){
s.diff = s.norm(coord=x,...)
l.s = lower.tri(s.diff,diag=FALSE)
sl = s.diff[l.s]
el = (1/(16*pi))*sl*log(sl)
d.s = !(l.s + t(l.s))
Em = matrix(0,nrow=nrow(s.diff),ncol=ncol(s.diff))
Em[l.s] = el
Em[t(l.s)] = el
Em
}
.tpsBase = function(x,tol=.Machine$double.neg.eps,...){
eta.r <- function(r) ifelse(r>0,(1/(16*pi))*r**2*log(r**2),0)
dist.titj <- .h2D(as.matrix(x),...)$h
aux.Q <- apply(dist.titj,2,eta.r)
list(V=diag(1,nrow(aux.Q)), A=solve(aux.Q,tol=tol), Ai=aux.Q, Asr=diag(1,nrow(aux.Q)), K=aux.Q, type="tps")
}
.csBase = function (x){
knots.values = .getKnots(x)
nv = length(x)
kv = length(knots.values)
hv = knots.values[2:kv] - knots.values[1:(kv - 1)]
Vm = matrix(0, nrow = kv, ncol = kv - 2)
Am = matrix(0, nrow = kv - 2, ncol = kv - 2)
for (k in 1:(kv - 2)) {
Vm[k, k] = 1/(hv[k])
Vm[k + 1, k] = -(1/(hv[k]) + 1/(hv[k + 1]))
Vm[k + 2, k] = 1/(hv[k + 1])
}
for (k in 1:(kv - 2)) {
if (k < kv - 2) {
Am[k, k + 1] = hv[k + 1]/6
Am[k + 1, k] = hv[k + 1]/6
Am[k, k] = (hv[k] + hv[k + 1])/3
}
if (k == (kv - 2)) {
Am[k, k] = (hv[k] + hv[k + 1])/3
}
}
KKm = Vm %*% qr.solve(Am) %*% t(Vm)
A.im.sr = svd(Am)
A.im.sr = A.im.sr$v %*% diag(sqrt(A.im.sr$d)) %*% t(A.im.sr$v)
list(V = Vm, A = Am, K = KKm, Asr = A.im.sr)
}
.psBase = function(x, diff.order = 1:1000000){
.defaultArg("diff.order",2)
x.tau = .getKnots(x)
aux.D = diff(diag(length(x.tau)), differences = diff.order)
aux.Q = t(aux.D)%*%aux.D
list(V=t(aux.D), A=diag(1,nrow(aux.D)), Asr=diag(1,nrow(aux.D)), K=aux.Q)
}
.tensorBase = function(obj=NULL, grid.list=NULL, type=c("cs","ps"), diff.order=2){
.defaultArg("type","cs")
if(is.null(obj) && is.null(grid.list)) stop("obj or grid.list must be specified!")
if(!is.null(grid.list) && !is.matrix(grid.list) && is.null(obj)) grid.list = grid.list
if(!is.null(grid.list) && is.matrix(grid.list) && is.null(obj)) grid.list = attributes(grid.list)$grid.list
if(!is.null(obj) && is.sss(obj)) grid.list = attributes(obj@grid)$grid.list
if(is.null(grid.list) && !is.null(obj) && is.sss(obj)) grid.list = obj@knots
II.x = diag(1,length(grid.list[[1]]))
JJ.y = diag(1,length(grid.list[[2]]))
if(type=="cs"){
QQ.x = .csBase(grid.list[[1]])
QQ.y = .csBase(grid.list[[2]])
} else {
QQ.x = .psBase(grid.list[[1]],diff.order)
QQ.y = .psBase(grid.list[[2]],diff.order)
}
list(Ix=II.x,Qx=QQ.x,Jy=JJ.y,Qy=QQ.y,type=type,difference=if(type=="ps"){diff.order})
}
.tensorPenalty = function(obj, alpha = c(1,1), type=c("Q","P","both")){
.defaultArg("type","P")
if(missing(obj) || is.null(obj)) stop("tensor base obj must be specified!")
if(length(alpha)!=2L || !is.numeric(alpha)) stop("alpha must be a 2*1 vector of smoothing parameters for the coordinates")
if(type=="Q") return(obj$Jy%x%obj$Qx$K + obj$Qy$K%x%obj$Ix)
if(type=="P") return(alpha[1]*(obj$Jy%x%obj$Qx$K) + alpha[2]*(obj$Qy$K%x%obj$Ix))
if(type=="both") list(Qxy = obj$Jy%x%obj$Qx$K + obj$Qy$K%x%obj$Ix, Pxy = alpha[1]*(obj$Jy%x%obj$Qx$K) + alpha[2]*(obj$Qy$K%x%obj$Ix))
}
.decomposePenalty = function(PM, method.eigen = TRUE, diagA = FALSE, thr = .Machine$double.eps, ftr = ifelse(method.eigen,1,10), inverse.method = c("solve","ginv",".invEigen",".invSchur",".invGsvd"), itol = .Machine$double.neg.eps){
.defaultArg("inverse.method","solve")
if(missing(PM)) stop("PM must be a penalty matrix!")
PM = as.matrix(PM)
if(nrow(PM)!=ncol(PM)) stop("PM must be a squared matrix!")
d.PM = do.call(ifelse(method.eigen,"eigen","svd"),list(x=PM))
d.rule = d.PM[[ifelse(method.eigen,"values","d")]] > sqrt(thr)*max(d.PM[[ifelse(method.eigen,"values","d")]])*ftr
if(!diagA){
V.el = d.PM[[ifelse(method.eigen,"vectors","u")]][,d.rule]
Ai.el = diag(d.PM[[ifelse(method.eigen,"values","d")]][d.rule])
A.el = .whichInverse(inverse.method, Ai.el, tol = itol, Re.values = TRUE)
A.sr.el = sqrt(A.el)
} else {
V.el = (d.PM[[ifelse(method.eigen,"vectors","u")]][,d.rule])%*%sqrt(diag(d.PM[[ifelse(method.eigen,"values","d")]][d.rule]))
Ai.el = diag(1,sum(d.rule))
A.el = Ai.el
A.sr.el = Ai.el
}
list(V = V.el, A = A.el, Ai = Ai.el, Asr = A.sr.el)
}
.generateXZ <- function(coords, QList){
if(missing(QList)) QList <- NULL
if(is.null(QList)) stop("QList must be specified!")
X1 <- cbind(1,.getKnots(coords[,1]))
X2 <- cbind(1,.getKnots(coords[,2]))
Z1 <- QList$Qx$V%*%solve(t(QList$Qx$V)%*%QList$Qx$V)
Z2 <- QList$Qy$V%*%solve(t(QList$Qy$V)%*%QList$Qy$V)
X <- cbind(X2%x%X1)
Z <- cbind(X2%x%Z1,Z2%x%X1,Z2%x%Z1)
list(X = X, Z = Z)
}
.generateV <- function(Z, itol = .Machine$double.neg.eps) return(Z%*%solve(t(Z)%*%Z, tol = itol))
.generateAi <- function(Z, Q, itol = .Machine$double.neg.eps) return(t(Z)%*%Q%*%Z)
.generateZs <- function(Rs, V, itol = .Machine$double.neg.eps) return(Rs%*%V%*%solve(t(V)%*%Rs%*%V, tol = itol))
.incidenceMatrix = function(x, by = NULL){
knots.values = .getKnots(x)
if(is.null(by)){
nv = length(x)
} else {
knots.i = by(x,by,.getKnots,simplify=TRUE)
nprod = 1
for(j in 1:length(by)){
nprod = nprod*nlevels(factor(by[[j]]))
}
nv = nprod*length(knots.values)
NAknots = rep(NA,length(knots.values))
NAn = NULL
for(j in 1:length(knots.i)){
NAn[[j]] = NAknots
NAn[[j]][match(knots.i[[j]],knots.values)] = knots.i[[j]]
}
NAn = unlist(NAn)
x = NAn
}
kv = length(knots.values)
W = matrix(0, nrow = nv, ncol = kv)
xc = paste(round(x, options()$digits))
kc = paste(round(knots.values, options()$digits))
for (i in 1:nv) {
for (r in 1:kv) {
if (xc[i] == kc[r]) {
W[i, r] = 1
}
}
}
return(W)
}
.incidenceMatrix = function(x, by = NULL, add = FALSE){
knots.values = .getKnots(x)
if(is.null(by)){
nv = length(x)
} else {
knots.i = by(x,by,.getKnots,simplify=TRUE)
nprod = 1
for(j in 1:length(by)){
nprod = nprod*nlevels(factor(by[[j]]))
}
nv = nprod*length(knots.values)
NAknots = rep(NA,length(knots.values))
NAn = NULL
for(j in 1:length(knots.i)){
NAn[[j]] = NAknots
NAn[[j]][match(knots.i[[j]],knots.values)] = knots.i[[j]]
}
NAn = unlist(NAn)
x = NAn
}
kv = length(knots.values)
W = matrix(0, nrow = nv, ncol = kv - ifelse(add, 1, 0))
xc = paste(round(x, options()$digits))
kc = paste(round(knots.values, options()$digits))
for (i in 1:nv) {
for (r in 1:(kv - ifelse(add, 1, 0))) {
if(add){
if (xc[i] == kc[r] && xc[i] != kc[kv]) {
W[i, r] = 1
} else if (xc[i] == kc[kv]) {
W[i, r] = -1
}
} else {
if (xc[i] == kc[r]) {
W[i, r] = 1
}
}
}
}
return(W)
}
.incidenceSpatial <- function(coords,grid,obs){
if(missing(grid)) grid <- NULL
if(missing(obs)) obs <- NULL
if(is.null(coords)) stop(".incidenceSpatial: coordinates must be defined!")
coords.id = apply(coords,1,paste,sep="",collapse=":")
if(is.null(grid)) grid <- .createGrid(coords)
grid.id = apply(grid,1,paste,sep="",collapse=":")
if(is.null(grid.knots <- attributes(grid)$grid.list)){
n.x = length(.getKnots(coords[,1]))
n.y = length(.getKnots(coords[,2]))
} else {
n.x = length(grid.knots[[1]])
n.y = length(grid.knots[[2]])
}
if(is.null(obs)) obs <- rep(1,length(coords.id))
row.id = match(coords.id,grid.id)
W.coords = matrix(0, nrow = n.x*n.y, ncol = length(coords.id))
for(i in 1:length(row.id)){
W.coords[row.id[i],i] = ifelse(!is.na(obs[i]),1,0)
}
W.coords <- t(W.coords)
attr.Wxy = attributes(W.coords)
attributes(W.coords) = attr.Wxy
return(W.coords)
}
.incidenceSpatial <- function(coords,grid,obs,add=FALSE){
if(missing(grid)) grid <- NULL
if(missing(obs)) obs <- NULL
if(is.null(coords)) stop(".incidenceSpatial: coordinates must be defined!")
coords.id = apply(coords,1,paste,sep="",collapse=":")
if(is.null(grid)) grid <- .createGrid(coords)
grid.id = apply(grid,1,paste,sep="",collapse=":")
if(is.null(grid.knots <- attributes(grid)$grid.list)){
n.x = length(.getKnots(coords[,1]))
n.y = length(.getKnots(coords[,2]))
} else {
n.x = length(grid.knots[[1]])
n.y = length(grid.knots[[2]])
}
if(is.null(obs)) obs <- rep(1,length(coords.id))
row.id = match(coords.id,grid.id)
W.coords = matrix(0, nrow = n.x*n.y, ncol = length(coords.id))
for(i in 1:length(row.id)){
W.coords[row.id[i],i] = ifelse(!is.na(obs[i]),1,0)
}
W.coords <- t(W.coords)
if(add){
W.i <- cbind(.incidenceMatrix(1:length(coords.id),add=TRUE),0)
W.coords = W.i%*%W.coords
}
attr.Wxy = attributes(W.coords)
attributes(W.coords) = attr.Wxy
return(W.coords)
}
.genZ = function(V, R = NULL, Asr = NULL, Ztr = c("V","VA","RV","RVA"), only.z = FALSE,
inverse.method = c("solve","ginv",".invEigen",".invSchur",".invGsvd"),
itol = .Machine$double.neg.eps){
.defaultArg("inverse.method","solve")
if(missing(V)) stop("Matrix V must be specified!")
.defaultArg("Ztr","RV")
if(Ztr=="RV" | Ztr=="RVA"){
VRV = t(V)%*%R%*%V
VRVi = .whichInverse(inverse.method, VRV, tol=itol, Re.values=TRUE)
Z = R%*%V%*%VRVi
}
if(Ztr=="RVA") Z = Z%*%Asr
if(Ztr=="V" | Ztr=="VA"){
VRV = t(V)%*%V
VRVi = .whichInverse(inverse.method, VRV, tol=itol, Re.values=TRUE)
Z = V%*%VRVi
}
if(Ztr=="VA") Z = Z%*%Asr
if(only.z) return(Z) else list(Z=Z,VRV=VRV,VRVi=VRVi)
}
.autoCor <- function(model = c("identity","pure.nugget","matern","powered.exponential","spherical","wave","exponential","gaussian","cubic","circular","gencauchy","cauchy","RMmatern","RMwhittle","RMgneiting","RMgengneiting","RMnugget","RMcauchy","RMexp","RMgencauchy","RMgauss","RMspheric","RMstable","RMpoweredexp"), what = c("expr","bounds","all"), ...){
m.c <- match.call(expand.dots=TRUE)
.defaultArg("model","RMwhittle")
.defaultArg("what","expr")
model.RF = c("RMmatern","RMwhittle","RMgneiting","RMgengneiting","RMnugget","RMcauchy","RMexp","RMgencauchy","RMgauss","RMspheric","RMstable","RMpoweredexp")
el = list()
if(what == "expr" | what == "all"){
add.args = list(...)
if(is.null(add.args$tol.nugget)) tol.nugget <- 1.0e-15 else tol.nugget <- add.args$tol.nugget
if(model == "matern") el$rho.h = expression(ifelse(h==0,1,((2**(1-kappa))/gamma(kappa))*((h/phi)**kappa)*besselK(h/phi, nu = kappa, expon.scaled = FALSE)))
if(model == "powered.exponential") el$rho.h = expression(exp(-(h/phi)**kappa))
if(model == "spherical") el$rho.h = expression(ifelse(h >= 0 & h <= phi,1-(3/2)*(h/phi)+(1/2)*((h/phi)**3),0))
if(model == "wave") el$rho.h = expression(ifelse(h>0,((h/phi)**(-1))*sin(h/phi),1))
if(model == "cauchy") el$rho.h = expression(ifelse(h>0,(1+(h/phi)**2)**(-kappa),1))
if(model == "gencauchy") el$rho.h = expression(ifelse(h>0,(1+(h/phi)**kappa[2])**(-kappa[1]/kappa[2]),1))
if(model == "circular") el$rho.h = expression(ifelse(h >= 0 & h < phi,1-2*(min(h/phi,1)*sqrt(1-min(h/phi,1)**2)+asin(sqrt(h/phi)))/pi,0))
if(model == "cubic") el$rho.h = expression(ifelse(h >= 0 & h < phi,1-(7*(h/phi)**2-8.75*(h/phi)**3+3.5*(h/phi)**5-0.75*(h/phi)**7),0))
if(model == "gaussian") el$rho.h = expression(exp(-(h/phi)**2))
if(model == "exponential") el$rho.h = expression(exp(-(h/phi)))
if(model == "identity") el$rho.h = expression(ifelse(h==0,1,0))
if(model == "pure.nugget") el$rho.h = expression(ifelse(h>=0 & h<tol.nugget,1,0))
if(!is.na(match(model,model.RF))){
requireNamespace("RandomFields", quietly = TRUE)
l.args = formals(get(model))
add.args = add.args[!is.na(match(names(add.args),names(l.args)))]
if(any(grepl("Aniso",names(add.args)))) add.args[[grepl("Aniso",names(add.args))]] <- m.c[["Aniso"]]
if(any(grepl("proj",names(add.args)))) add.args[[grepl("proj",names(add.args))]] <- m.c[["proj"]]
plug.r = ifelse(length(add.args)>0,paste(",",paste(names(add.args),add.args,sep="=",collapse=","),sep=""),"")
plug.f = ifelse(length(add.args)>0,paste(",",paste(names(add.args),names(add.args),sep="=",collapse=","),sep=""),"")
r.args = paste("phi",ifelse(any(!is.na(match(model,c("RMmatern","RMwhittle","RMgengneiting","RMcauchy","RMgencauchy","RMstable","RMpoweredexp")))),",kappa",""),plug.r,sep="")
f.args = paste("var=1,scale=phi",
ifelse(any(!is.na(match(model,c("RMmatern","RMwhittle")))),",nu=kappa",
ifelse(any(!is.na(match(model,c("RMgengneiting")))),",mu=kappa[1],kappa=kappa[2]",
ifelse(any(!is.na(match(model,c("RMcauchy")))),",gamma=kappa",
ifelse(any(!is.na(match(model,c("RMgencauchy")))),",beta=kappa[1],alpha=kappa[2]",
ifelse(any(!is.na(match(model,c("RMstable","RMpoweredexp")))),",alpha=kappa",""))))),
plug.f,sep="")
}
if(any(!is.na(match(model,model.RF)))){
el$rho.h = eval(parse(text=paste("function(",r.args,"){",paste("RandomFields::",model,sep=""),"(",f.args,")}",sep="")))
environment(el$rho.h) <- parent.frame()
}
}
if(what == "bounds" | what == "all"){
if(!is.na(match(model,c("matern","cauchy","RMmatern","RMwhittle","RMcauchy")))){
el$lower = c(phi=0,kappa=0)
el$upper = c(phi=Inf,kappa=Inf)
}
if(!is.na(match(model,c("RMgengneiting")))){
el$lower = c(phi=0,kappa1=1,kappa2=0)
el$upper = c(phi=Inf,kappa1=Inf,kappa2=3)
}
if(!is.na(match(model,c("gencauchy","RMgencauchy")))){
el$lower = c(phi=0,kappa1=0,kappa2=0)
el$upper = c(phi=Inf,kappa1=Inf,kappa2=2)
}
if(!is.na(match(model,c("powered.exponential","RMstable","RMpoweredexp")))){
el$lower = c(phi=0,kappa=0)
el$upper = c(phi=Inf,kappa=2)
}
if(!is.na(match(model,c("identity","pure.nugget","spherical","wave","exponential","gaussian","cubic","circular","RMnugget","RMgneiting","RMexp","RMgauss","RMspheric")))){#"RMnugget"?
el$lower = c(phi=0)
el$upper = c(phi=Inf)
}
}
if(what == "expr"){
return(el$rho.h)
} else {
return(el)
}
}
.semiVar = function(model = c("matern","gaussian","exponential","power","cubic","penta.spherical","spherical","wave","sin.hole","identity","pure.nugget"), what = c("expr","bounds","all"), tol.nugget = 1.0e-15){
.defaultArg("model","matern")
.defaultArg("what","expr")
sl = list()
if(what == "expr" | what == "all"){
if(model == "matern") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill*(1-(2**(1-kappa))/gamma(kappa))*((h/phi)**kappa)*besselK(h/phi, nu = kappa, expon.scaled = FALSE),nugget+0))
if(model == "gaussian") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill*(1-exp(-(h/phi)**2)),nugget+0))
if(model == "exponential") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill*(1-exp(-(h/phi))),nugget+0))
if(model == "power") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill*(h**phi),nugget+0))
if(model == "cubic") sl$gamma.h = expression(ifelse((h > tol.nugget & h <= phi),nugget+sill*(7*((h/phi)**2)-(35/4)*((h/phi)**3)+(7/2)*((h/phi)**5)-(3/4)*((h/phi)**7)),ifelse(tol.nugget<phi & phi<h,nugget+sill,nugget+0)))
if(model == "penta.spherical") sl$gamma.h = expression(ifelse((h > tol.nugget & h <= phi),nugget+sill*((15/8)*(h/phi)-(5/4)*((h/phi)**3)+(3/8)*((h/phi)**5)),ifelse(tol.nugget<phi & phi < h,nugget+sill,nugget+0)))
if(model == "spherical") sl$gamma.h = expression(ifelse((h > tol.nugget & h <= phi),nugget+sill*((3/2)*(h/phi)+(1/2)*((h/phi)**3)),ifelse(tol.nugget<phi & phi < h,nugget+sill,nugget+0)))
if(model == "wave") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill*(1-((h/phi)**(-1))*sin(h/phi)),nugget+0))
if(model == "sin.hole") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill*(1-((pi*h/phi)**(-1))*sin(pi*h/phi)),nugget+0))
if(model == "pure.nugget") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+0,0))
if(model == "identity") sl$gamma.h = expression(ifelse(h > tol.nugget,nugget+sill,nugget+0))
}
if(what == "bounds" | what == "all"){
if(model == "matern"){
sl$lower = c(sill=0, phi=0, kappa=0, nugget=0)
sl$upper = c(sill=Inf, phi=Inf, kappa=Inf, nugget=Inf)
}
if(!is.na(match(model, c("gaussian","exponential","power","cubic","penta.spherical","spherical","wave","sin.hole")))){
sl$lower = c(sill=0, phi=0, nugget=0)
sl$upper = c(sill=Inf, phi=Inf, nugget=Inf)
}
if(model == "identity"){
sl$lower = c(sill=0, nugget=0)
sl$upper = c(sill=Inf, nugget=Inf)
}
if(model == "pure.nugget"){
sl$lower = c(nugget=0)
sl$upper = c(nugget=Inf)
}
}
if(what == "expr"){
return(sl$gamma.h)
} else {
return(sl)
}
}
.autoCorF = function(h, phi = 1, kappa = NULL, model = c("matern","identity","pure.nugget","powered.exponential","spherical","wave","exponential","gaussian","cubic","circular","gencauchy","cauchy","RMmatern","RMwhittle","RMgneiting","RMgengneiting","RMnugget","RMcauchy","RMexp","RMgencauchy","RMgauss","RMspheric","RMstable","RMpoweredexp"), ...){
m.c <- match.call(expand.dots=TRUE)
if(missing(phi) | is.null(phi)){phi = 1}
.defaultArg("model","RMwhittle")
add.args = list(...)
if(is.null(add.args$tol.nugget)) tol.nugget <- 1.0e-15 else tol.nugget <- add.args$tol.nugget
model.RF = c("RMmatern","RMwhittle","RMgneiting","RMgengneiting","RMnugget","RMcauchy","RMexp","RMgencauchy","RMgauss","RMspheric","RMstable","RMpoweredexp")
if(!is.na(match(model,model.RF))){
requireNamespace("RandomFields", quietly = TRUE)
}
rho = .autoCor(model = model, what = "all", ...)
if(any(grepl("kappa",names(rho$lower))) & is.null(kappa)){
stop("Parameter kappa was not defined!")
}
theta = c(phi=phi,kappa=kappa)
for(i in 1:length(rho$lower)){
if(!(theta[i]>=rho$lower[i] && theta[i]<=rho$upper[i])){
stop(paste("Parameter ",names(theta)[i]," lies outside allowed region for its family!\nEnter ",names(theta)[i]," value in ",ifelse(rho$lower[i]==0,"(","["),rho$lower[i],",",rho$upper[i],ifelse(rho$upper[i]==Inf,").","]."),sep="",collapse=""))
}
}
if(!any(grepl("kappa",names(rho$lower))) & !is.null(kappa)){
message("kappa not needed for the family, so not used!")
}
h.mat = (is.matrix(h) || is.data.frame(h)) && nrow(h)==ncol(h)
if(h.mat){
if(!isSymmetric(h) || nrow(h)!=ncol(h)) stop("h must be a vector or a symmetric n*n matrix of distances!")
Rho.h = matrix(1,nrow=nrow(h),ncol=ncol(h))
low.h = lower.tri(Rho.h,diag=FALSE)
h = h[low.h]
}
if(any(grepl(model,c("matern","identity","pure.nugget","powered.exponential","spherical","wave","exponential","gaussian","cubic","circular","gencauchy","cauchy")))){
aux = eval(rho$rho.h)
} else {
call.cf = eval(rho$rho.h)
fargs <- formals(call.cf)
nargs <- names(formals(call.cf))
largs <- list()
for(i in 1:length(nargs)){
if(any(grepl(nargs[i],names(theta)))){
largs[[i]] <- theta[grepl(nargs[i],names(theta))]
} else {
largs[[i]] <- fargs[[i]]
}
names(largs)[i] <- nargs[i]
}
d.c <- call("call.cf")
for(i in 1:length(largs)){
if(names(largs)[i]=="Aniso"){
d.c[[names(largs)[i]]] <- m.c[["Aniso"]]
} else if(names(largs)[i]=="proj") {
d.c[[names(largs)[i]]] <- m.c[["proj"]]
} else {
d.c[[names(largs)[i]]] <- largs[[i]]
}
}
e.c <- eval(d.c)
environment(e.c) <- parent.frame()
aux = RandomFields::RFcov(e.c, distances=h, dim=2)
}
if(h.mat){
Rho.h[low.h] = aux
Rho.h = t(Rho.h)
Rho.h[low.h] = aux
return(Rho.h)
} else {
return(aux)
}
}
.semiVarF = function(h, phi = NULL, sill = NULL, kappa = NULL, nugget = NULL, model = c("matern","gaussian","exponential","power","powered.exponential","cubic","circular","penta.spherical","spherical","wave","sin.hole","cauchy","gencauchy","identity","pure.nugget","RMmatern","RMwhittle","RMgneiting","RMgengneiting","RMnugget","RMcauchy","RMexp","RMgencauchy","RMgauss","RMspheric","RMstable","RMpoweredexp"), use.cor = FALSE, tol.nugget = 1.0e-15, ...){
m.c = match.call(expand.dots = TRUE)
cor.model = paste(formals(.autoCorF)$model)[-1]
only.var = paste(formals(.semiVar)$model)[-1]
if(missing(use.cor)) use.cor <- NULL
if(is.null(use.cor)){
in.cor <- any(grepl(model,cor.model))
in.var <- any(grepl(model,only.var))
use.cor <- ifelse((in.cor & in.var) | (in.cor & !in.var), TRUE,ifelse(!in.cor & in.var, FALSE, NA))
}
if(missing(phi) | is.null(phi)){phi = NULL}
if(missing(sill) | is.null(sill)){sill = NULL}
if(missing(kappa) | is.null(kappa)){kappa = NULL}
if(missing(nugget) | is.null(nugget)){nugget = 0}
model = match.arg(model)
model.RF = c("RMmatern","RMwhittle","RMgneiting","RMgengneiting","RMnugget","RMcauchy","RMexp","RMgencauchy","RMgauss","RMspheric","RMstable","RMpoweredexp")
if(!is.na(match(model,model.RF))){
requireNamespace("RandomFields", quietly = TRUE)
}
if(use.cor){
rho.v = .autoCorF(h = h, phi = phi, kappa = kappa, model = model, ...)
return(as.vector(nugget)*(1-(h<tol.nugget)) + as.vector(sill)*(1 - rho.v))
} else {
gamma.h = .semiVar(model = model, what = "all", tol.nugget = tol.nugget)
if(any(grepl("kappa",names(gamma.h$lower))) & is.null(kappa)){
stop("Parameter kappa was not defined!")
}
theta.all = c(sill=as.vector(sill),phi=as.vector(phi),kappa=as.vector(kappa),nugget=as.vector(nugget))
isnt = is.na(match(names(theta.all),names(gamma.h$lower)))
theta = theta.all[!isnt]
not.given = is.na(match(names(gamma.h$lower),names(theta)))
if(any(not.given)){
stop(paste("Parameter",paste(names(gamma.h$lower)[not.given],sep="",collapse=", "),"not given but needed for the covariance family!"))
}
if(any(isnt)){
message(paste("Parameter",paste(names(theta.all)[isnt],sep="",collapse=", "),"not needed for the family, so not used!"))
}
for(i in 1:length(gamma.h$lower)){
if(!(theta[i]>=gamma.h$lower[i] && theta[i]<=gamma.h$upper[i])){
stop(paste("Parameter ",names(theta)[i]," lies outside allowed region for its family!\nEnter ",names(theta)[i]," value in ",ifelse(gamma.h$lower[i]==0,"(","["),gamma.h$lower[i],",",gamma.h$upper[i],ifelse(gamma.h$upper[i]==Inf,").","]."),sep="",collapse=""))
}
}
h.mat = (is.matrix(h) || is.data.frame(h)) && nrow(h)==ncol(h)
if(h.mat){
if(!isSymmetric(h) || nrow(h)!=ncol(h)) stop("h must be a vector or a symmetric n*n matrix of distances!")
G.h = matrix(nugget,nrow=nrow(h),ncol=ncol(h))
low.h = lower.tri(G.h,diag=FALSE)
h = h[low.h]
}
if(any(grepl(model,c("matern","gaussian","exponential","power","cubic","penta.spherical","spherical","wave","sin.hole","identity","pure.nugget")))){
aux = eval(gamma.h$gamma.h)
} else {
call.cf = eval(gamma.h$gamma.h)
fargs <- formals(call.cf)
nargs <- names(formals(call.cf))
largs <- list()
for(i in 1:length(nargs)){
if(any(grepl(nargs[i],names(theta)))){
largs[[i]] <- theta[grepl(nargs[i],names(theta))]
} else {
largs[[i]] <- fargs[[i]]
}
names(largs)[i] <- nargs[i]
}
d.c <- call("call.cf")
for(i in 1:length(largs)){
if(names(largs)[i]=="Aniso"){
d.c[[names(largs)[i]]] <- m.c[["Aniso"]]
} else if(names(largs)[i]=="proj") {
d.c[[names(largs)[i]]] <- m.c[["proj"]]
} else {
d.c[[names(largs)[i]]] <- largs[[i]]
}
}
e.c <- eval(d.c)
environment(e.c) <- parent.frame()
aux = RandomFields::RFvariogram(e.c, distances=h, dim=2)
}
if(h.mat){
G.h[low.h] = aux
G.h = t(G.h)
G.h[low.h] = aux
return(G.h)
} else {
return(aux)
}
}
}
.h2D <- function(x, angle = 0, ratio = 1){
if(angle != 0 | ratio != 1){
H.m = matrix(c(cos(angle),sin(angle),-sin(angle),cos(angle)),2,2,byrow = TRUE)
D.m = diag(c(1,1/ratio))
A.m = D.m%*%H.m
} else {
A.m = diag(1,2)
}
.hNormRow <- function(coord,r.loc, A = NULL){
d.ro = NULL
for(r.other in 1:nrow(coord)){
h.diff = coord[r.loc,] - coord[r.other,]
if(!is.null(A)){
h.diff = A%*%h.diff
}
d.ro[r.other] = sqrt(t(h.diff)%*%h.diff)
}
return(as.numeric(d.ro))
}
hs.row = NULL
for(r.row in 1:nrow(x)){
hs.row = cbind(hs.row,.hNormRow(x,r.row,A=A.m))
}
list(coord = x, h = hs.row)
}
.h2D <- function(x,angle=0,ratio=1) {
if (!is.numeric(x))
stop("argument x must be numeric")
if (!(nrow(x)>1) | !(ncol(x)==2))
stop("argument x must be of dimension n*2")
if (!(ratio > 0 & ratio <= 1) | !(angle >= 0 & angle <= 360))
stop("ratio must belong to (0,1] and angle to [0,360]")
n = as.integer(nrow(x))
s = apply(x,2,as.double)
a = as.double(angle)
r = as.double(ratio)
h = matrix(as.double(rep(0,n*n)),n,n)
out <- .Fortran("h2d", n=n, s=s, a=a, r=r, h=h, PACKAGE = "scpm")
list(coord=x,h=out$h)
}
.createCov = function(dist.m = NULL, coords.m = NULL, to.Inf = 1.0e+24, rhos = NULL, ...){
force.cor <- FALSE
if(missing(rhos)) rhos <- NULL
if(missing(to.Inf)) to.Inf = 1.0e+24
if(is.null(to.Inf)) to.Inf = 1.0e+24
add.list = list(...)
if(!is.null(add.list$phi)){ phi = add.list$phi } else { phi = NULL }
if(!is.null(add.list$sill)){ sill = add.list$sill } else {sill = NULL }
if(!is.null(add.list$kappa)){ kappa = add.list$kappa } else { kappa = NULL }
if(!is.null(add.list$nugget)){ nugget = add.list$nugget } else { nugget = NULL }
if(!is.null(add.list$model)){ model = add.list$model } else { model = NULL }
if(is.null(add.list$use.cor)){ use.cor = FALSE } else { use.cor = add.list$use.cor }
if(is.null(add.list$tol.nugget)){ tol.nugget = 1.0e-15 } else { tol.nugget = add.list$tol.nugget }
if(is.null(nugget) & is.null(sill) & is.null(rhos)){
stop(".createCov: nugget,sill or rhos must be defined!")
} else if(!is.null(nugget) & !is.null(sill) & is.null(rhos)){
rhos <- sill/(nugget + sill)
} else if(!is.null(nugget) & is.null(sill) & !is.null(rhos)){
sill <- nugget/(1-rhos) - nugget
} else if(is.null(nugget) & !is.null(sill) & !is.null(rhos)){
nugget <- sill/rhos - sill
} else if(is.null(nugget) & is.null(sill) & !is.null(rhos)){
force.cor <- TRUE
} else {
stop(".createCov: nugget and sill, only rhos or rhos with nugget or sill must be specified!")
}
if(missing(dist.m)) dist.m <- NULL
if(missing(coords.m)) coords.m <- NULL
if(is.null(add.list$angle)) { angle = 0 } else { angle = add.list$angle }
if(is.null(add.list$ratio)) { ratio = 1 } else { ratio = add.list$ratio }
lcor = paste(formals(.autoCorF)$model)[-1]
lvar = paste(formals(.semiVarF)$model)[-1]
if(is.null(dist.m) & is.null(coords.m)){
stop("matrix of distances or coordinates must be given!")
} else if(is.null(dist.m) & !is.null(coords.m)){
dist.m <- .h2D(as.matrix(coords.m), angle = angle, ratio = ratio)$h
}
in.cor = any(!is.na(match(add.list$model,lcor)))
in.var = any(!is.na(match(add.list$model,lvar)))
use.var = ifelse(((in.cor & in.var) | (!in.cor & in.var)) & !force.cor,TRUE,ifelse((in.cor & in.var) & force.cor,FALSE,ifelse((!in.cor & in.var) & force.cor,TRUE,ifelse(in.cor & !in.var,FALSE,NA))))
if(!is.na(use.var)){
if(use.var){
if(ifelse(in.var & in.cor & !use.cor,TRUE,FALSE)){
use.cor = !use.cor
}
if(ifelse(in.var & !in.cor & use.cor,TRUE,FALSE)){
use.cor = !use.cor
}
if(force.cor) stop(".createCov: model does not allow only rhos. Additional nugget or sill must be specified!")
c.gI = call(".semiVarF",h=to.Inf,phi=phi,sill=sill,kappa=kappa,nugget=0,model=model,use.cor=use.cor)
gamma.Inf = eval(c.gI)
c.gh = call(".semiVarF",h=dist.m,phi=phi,sill=sill,kappa=kappa,nugget=0,model=model,use.cor=use.cor)
gamma.h = eval(c.gh)
c.mat = gamma.Inf - gamma.h
gamma.n = call(".semiVarF",h=dist.m,phi=phi,sill=0,kappa=kappa,nugget=nugget,model="pure.nugget",tol.nugget=tol.nugget,use.cor=FALSE)
c.mat = eval(gamma.n) + c.mat
c.mat[c.mat==Inf] = to.Inf
c.mat[c.mat==-Inf] = -to.Inf
return(c.mat)
} else {
if(!force.cor){
sills <- nugget + sill
c.gh = call(".semiVarF",h=dist.m,nugget=nugget,model="pure.nugget",tol.nugget=tol.nugget,use.cor=FALSE)
c.rh = call(".autoCorF",h=dist.m,phi=phi,kappa=kappa,model=model)
c.mat = (1/sills)%x%eval(c.gh) + rhos%x%eval(c.rh)
} else {
c.rh = call(".autoCorF",h=dist.m,phi=phi,kappa=kappa,model=model)
c.mat <- (1-rhos)%x%(dist.m<tol.nugget) + rhos%x%eval(c.rh)
}
c.mat[c.mat==Inf] = to.Inf
c.mat[c.mat==-Inf] = -to.Inf
return(c.mat)
}
} else {
stop("model not in the list of available models!")
}
}
.scpLoglik <- function(par, beta.par = FALSE, is.log.rhos = TRUE,
obs.z, dist.m, W, Xs, Z, V, fix.nugget = NULL, fix.kappa = NULL,
pty = "none", Q.base = NULL, Ztransf = "RV", model, ml = TRUE, use.profile = TRUE,
inverse.method = "solve", itol = .Machine$double.eps, createA = "eigen",
diagA = TRUE, nugget.tol = 1.0e-15, print.pars = TRUE, ...){
ac <- paste(formals(.autoCor)$model)[-1]
sv <- paste(formals(.semiVar)$model)[-1]
ul <- unique(c(ac,sv))
fl <- sv[is.na(match(sv,ac))]
non.stat <- any(grepl(paste("^",model,"$",sep="",collapse=""),fl))
nugget.par <- ifelse(non.stat & is.null(fix.nugget),TRUE,FALSE)
fit.nugget = ifelse(is.null(fix.nugget),TRUE,FALSE)
fit.rhos = ifelse(fit.nugget,TRUE,FALSE)
fit.kappa = ifelse(is.null(fix.kappa),TRUE,FALSE)
drop.nugget = ifelse(!non.stat | (non.stat & !fit.nugget),TRUE,FALSE)
drop.rhos = ifelse((!non.stat & !fit.nugget) | (non.stat & !fit.nugget),TRUE,FALSE)
n.xy = length(obs.z)
if(!beta.par) add.beta = 0 else add.beta = ncol(Xs)
if(beta.par) beta.s = par[1:add.beta] else beta.s <- NULL
if(pty=="none") smooth <- FALSE else smooth <- TRUE
if(smooth){
if(pty=="cs" | pty=="ps"){
add.alpha <- 2
} else if(pty=="tps"){
add.alpha <- 1
} else if(pty=="none"){
add.alpha <- 0
} else { stop("Non valid smooth type definition!") }
} else {
add.alpha <- 0
}
if(smooth){
alpha.xy = par[(add.beta+1):(add.beta+add.alpha)]
alpha <- sum(alpha.xy)
rho.xy = alpha.xy/alpha
} else {
alpha.xy = NULL
alpha <- NULL
rho.xy = NULL
}
if(!drop.rhos){
next.pos = add.beta + add.alpha + 2
if(is.log.rhos){
log.rhos <- c(log.rhos = par[add.beta + add.alpha + 1]) #log(sigma2/sigma2s)
rhos <- c(rhos = exp(log.rhos)) #sigma2/sigma2s
} else {
rhos <- c(rhos = par[add.beta + add.alpha + 1]) #sigma2/sigma2s
}
} else {
if(!is.null(fix.nugget)){
next.pos = add.beta + add.alpha + 1
rhos <- c(rhos = 1 - fix.nugget)
log.rhos <- c(log.rhos = log(rhos))
} else {
stop("A value for fix.nugget must be specified!")
}
}
isI = ifelse(model=="identity",TRUE,FALSE)
if(!isI){
phi = par[next.pos]
if(any(model==c("matern", "powered.exponential", "cauchy", "gencauchy", "RMmatern", "RMwhittle", "RMgengneiting", "RMcauchy", "RMgencauchy", "RMstable", "RMpoweredexp"))){
if(fit.kappa){ kappa = par[(next.pos+1):ifelse(fit.nugget & !drop.nugget,length(par)-1,length(par))] } else { kappa = fix.kappa }
} else {
kappa = NULL
}
} else {
phi = NULL
kappa = NULL
}
diff.length <- length(par)-length(c(beta.s,alpha.xy,if(fit.rhos & !drop.rhos){rhos},phi,if(is.null(fix.kappa)){kappa}))
if(sum((fit.nugget & !drop.nugget))!=sum(diff.length)){
if(sum((fit.nugget & !drop.nugget))>sum(diff.length)) stop(paste("Length of parameter defined different than required. Perhaps parameter tau2 must be included in par for model ",model,".",sep=""))
if(sum((fit.nugget & !drop.nugget))<sum(diff.length)) stop(paste("Length of parameter defined different than required. Perhaps parameter tau2 should not be included in par for model ",model,".",sep=""))
}
if(nugget.par){
tau2 <- c(tau2 = par[length(par)])
} else if(!is.null(fix.nugget)){
tau2 <- c(tau2 = fix.nugget)
} else {
tau2 <- 0
}
if(!non.stat){
R.s = .createCov(dist.m = dist.m, phi = phi, sill = NULL, kappa = kappa, nugget = NULL, rhos = rhos,
model = model, tol.nugget = nugget.tol, ...)
R.u = W%*%R.s%*%t(W)
} else {
sigma2s <- c(sigma2s = tau2/rhos)
sigma2 <- c(sigma2 = sigma2s - tau2)
Sigma.s = .createCov(dist.m = dist.m, phi = phi, sill = NULL, kappa = kappa, nugget = tau2, rhos = rhos,
model = model, tol.nugget = nugget.tol, ...)
R.s = (1/sigma2s)%x%Sigma.s
Sigma.u = W%*%Sigma.s%*%t(W)
R.u = (1/sigma2s)%x%Sigma.u
}
if(smooth){
if(pty=="tps"){
l.Q <- Q.base
Z <- Q.base$K
G.m = (1/alpha)%x%l.Q$A
Gi.m = alpha%x%l.Q$Ai
} else if(pty=="cs" | pty=="ps"){
Q.rho = .tensorPenalty(obj = Q.base, alpha = rho.xy, type = "P")
l.Q <- list(V = V)
if(Ztransf=="V" | Ztransf=="VA"){
Z <- Z
} else {
Z <- .generateZs(R.s, V, itol = itol)
}
l.Q$Ai <- .generateAi(Z, Q.rho, itol = itol)
l.Q$A <- solve(l.Q$Ai , tol = itol)
if(Ztransf=="VA" | Ztransf=="RVA"){
el.A <- eigen(l.Q$A)
l.Q$Asr <- el.A$vectors%*%diag(sqrt(el.A$values))%*%t(el.A$vectors)
Z <- Z%*%l.Q$Asr
G.m = (1/alpha)%x%diag(1,nrow(l.Q$A))
Gi.m = alpha%x%diag(1,nrow(l.Q$A))
} else {
G.m = (1/alpha)%x%l.Q$A
Gi.m = alpha%x%l.Q$Ai
}
}
Zs <- W%*%Z
LL = Zs%*%G.m%*%t(Zs) + R.u
} else {
LL = R.u
}
LL.i = .whichInverse(inverse.method, LL, tol = itol, Re.values = TRUE)
if(!beta.par | !ml) XLX = t(Xs)%*%LL.i%*%Xs
if(!beta.par){
XLXi = .whichInverse(inverse.method,XLX, tol = itol, Re.values = TRUE)
XLz = t(Xs)%*%LL.i%*%obs.z
beta.s = XLXi%*%XLz
}
e.c = (obs.z - Xs%*%beta.s)
RSS = t(e.c)%*%LL.i%*%e.c
den = ifelse(ml,n.xy,n.xy-ncol(Xs))
sigma2s = RSS/den
tau2 <- (1-rhos)*sigma2s
sigma2 = sigma2s - tau2
if(ml){
if(!use.profile){
logL = -0.5*den*log(2*pi) -0.5*determinant(sigma2s%x%LL)$modulus -0.5*RSS
}
plogL = -0.5*den*log(RSS) -0.5*determinant(LL)$modulus
} else {
if(!use.profile){
logL = -0.5*den*log(2*pi) -0.5*determinant(sigma2s%x%LL)$modulus -0.5*RSS -0.5*determinant(sigma2s%x%XLX)$modulus
}
plogL = -0.5*den*log(RSS) -0.5*determinant(LL)$modulus -0.5*determinant(XLX)$modulus
}
if(print.pars){
message("##########################################")
message("####### estimated parameter values #######")
message("##########################################")
message("##########################################")
if(use.profile){
message(paste("profile-logLik =",round(plogL,6)))
} else {
message(paste("logLik =",round(logL,6)))
}
message(paste(paste("beta.s[",1:length(beta.s),"] = ",sep=""),round(beta.s,6),sep="",collapse=" "))
if(exists("alpha.xy")) { if(!is.null(alpha.xy)) message(paste(paste("alpha[",1:length(alpha.xy),"] = ",sep=""),round(alpha.xy,6),sep="",collapse=" ")) }
if(exists("alpha")) { if(!is.null(alpha)) message(paste("alpha =",round(alpha,6))) }
if(exists("rhos")) { if(!is.null(rhos)) message(paste("rhos =",round(rhos,6))) }
if(exists("sigma2s")) { if(!is.null(sigma2s)) message(paste("sigma2s =",round(sigma2s,6))) }
if(exists("sigma2")) { if(!is.null(sigma2)) message(paste("sigma2 =",round(sigma2,6))) }
if(exists("phi")) { if(!is.null(phi)) message(paste("phi =",round(phi,6))) }
if(exists("kappa")){
if(class(kappa)!="function" & !is.null(kappa)){
message(paste(paste("kappa[",1:length(kappa),"] =",sep=""),round(kappa,6),sep="",collapse=" "))
}
}
if(exists("tau2")) { if(!is.null(tau2)) message(paste("tau2 =",round(tau2,6))) }
message("##########################################")
}
if(use.profile){
return(as.numeric(-plogL))
} else {
return(as.numeric(-logL))
}
}
.randomParsGrid = function(pars.bounds, control = list()){
if(is.null(control$is.zero)){ control$is.zero = 1.0e-6 }
if(is.null(control$is.inf)){ control$is.inf = 10000 }
if(is.null(control$size.each.grid)){ control$size.each.grid = 30 }
if(is.null(control$size.sample)){ control$size.sample = 25 }
c.zero = control$is.zero; c.Inf = control$is.inf; c.one = 1-control$is.zero
c.ninits = control$size.each.grid; c.lsize = control$size.sample
c.grid.0toInf = c(seq(c.zero,1,l=0.4*c.ninits), seq(1,10,l=0.4*c.ninits),
seq(10,100,l=0.15*c.ninits), seq(100,1000,l=0.025*c.ninits), seq(1000,c.Inf,l=0.025*c.ninits))
c.grid.0toInf = c.grid.0toInf[c.grid.0toInf>=c.zero & c.grid.0toInf<=c.Inf]
c.grid.mInfto0 = rev(-c.grid.0toInf)
c.grid.mInftoInf = c(seq(c.zero,1,l=0.4*c.ninits*0.5),seq(1,10,l=0.4*c.ninits*0.5),
seq(10,100,l=0.15*c.ninits*0.5),seq(100,1000,l=0.025*c.ninits*0.5),seq(1000,c.Inf,l=0.025*c.ninits*0.5))
c.grid.mInftoInf = c.grid.mInftoInf[c.grid.mInftoInf > c.zero & c.grid.mInftoInf < c.Inf]
c.grid.mInftoInf = c(-rev(c.grid.mInftoInf),c.grid.mInftoInf)
c.grid.0to2 = c(seq(c.zero,1,l=0.5*c.ninits),seq(1,2,l=0.5*c.ninits))
c.grid.0to1 = c(seq(c.zero,1,l=c.ninits))
c.pars.grid = list()
for(i in 1:length(pars.bounds$lower)){
if(pars.bounds$lower[i]==-Inf & pars.bounds$upper[i]==Inf) c.pars.grid[[i]] = unique(c.grid.mInftoInf)
if(pars.bounds$lower[i]==0 & pars.bounds$upper[i]==Inf) c.pars.grid[[i]] = unique(c.grid.0toInf)
if(pars.bounds$lower[i]==-Inf & pars.bounds$upper[i]==0) c.pars.grid[[i]] = unique(c.grid.mInfto0)
if(pars.bounds$lower[i]==0 & pars.bounds$upper[i]==2) c.pars.grid[[i]] = unique(c.grid.0to2)
if(pars.bounds$lower[i]==0 & pars.bounds$upper[i]==1) c.pars.grid[[i]] = unique(c.grid.0to1)
}
c.pars.grid = expand.grid(c.pars.grid)
colnames(c.pars.grid) = names(pars.bounds$lower)
c.ids = sample(1:nrow(c.pars.grid),min(nrow(c.pars.grid),c.lsize),replace=FALSE)
c.min = as.matrix(c.pars.grid[c.ids,])
colnames(c.min) = colnames(c.pars.grid)
rownames(c.min) = rownames(c.pars.grid)[c.ids]
return(c.min)
}
.getBounds = function(pars = list(fix.nugget = NULL, fix.kappa = NULL, is.log.rhos = TRUE, model = "matern"), smooth.type = c("none","cs","ps","tps")){
.defaultArg("smooth.type","none")
if(is.null(pars)) stop("pars list must be specified!")
if(is.null(pars$is.log.rhos)) pars$is.log.rhos = TRUE
ac <- paste(formals(.autoCor)$model)[-1]
sv <- paste(formals(.semiVar)$model)[-1]
ul <- unique(c(ac,sv))
fl <- sv[is.na(match(sv,ac))]
non.stat <- any(grepl(paste("^",pars$model,"$",sep="",collapse=""),fl))
nugget.par <- ifelse(non.stat & is.null(pars$fix.nugget),TRUE,FALSE)
fit.nugget = ifelse(is.null(pars$fix.nugget),TRUE,FALSE)
fit.rhos = ifelse(fit.nugget,TRUE,FALSE)
fit.kappa = ifelse(is.null(pars$fix.kappa),TRUE,FALSE)
c.bo.sv = .tryCatchWE(.semiVar(pars$model,"bounds"))
c.bo.ac = .tryCatchWE(.autoCor(pars$model,"bounds"))
c.use.var = ifelse((!c.bo.sv$error & !c.bo.ac$error) | (!c.bo.sv$error & c.bo.ac$error),
TRUE,
ifelse((c.bo.sv$error & !c.bo.ac$error),FALSE,NA))
if(!is.na(c.use.var)){
if(c.use.var){
c.lower = c.bo.sv$value$lower
c.upper = c.bo.sv$value$upper
} else {
c.lower = c.bo.ac$value$lower
c.upper = c.bo.ac$value$upper
}
} else {
stop("Model not in the list of models!")
}
c.names = names(c.lower)
if(any(grepl("sill",c.names))){
c.id <- grepl("sill",c.names)
c.names <- c.names[!c.id]
c.lower <- c.lower[!c.id]
c.upper <- c.upper[!c.id]
}
c.names = c(ifelse(pars$is.log.rhos,"log.rhos","rhos"),c.names)
c.lower = c(ifelse(pars$is.log.rhos,-Inf,0),c.lower)
c.upper = c(ifelse(pars$is.log.rhos,0,1),c.upper)
if(smooth.type=="cs" | smooth.type=="ps"){
c.names = c("alpha1","alpha2",c.names)
c.lower = c(0,0,c.lower)
c.upper = c(Inf,Inf,c.upper)
} else if(smooth.type=="tps"){
c.names = c("alpha",c.names)
c.lower = c(0,c.lower)
c.upper = c(Inf,c.upper)
}
drop.nugget = ifelse(!non.stat | (non.stat & !fit.nugget),TRUE,FALSE)
if(any(grepl("nugget",c.names))){
if(!drop.nugget){
c.names[grepl("nugget",c.names)] = "tau2"
} else {
c.id = grepl("nugget",c.names)
c.names = c.names[!c.id]
c.lower = c.lower[!c.id]
c.upper = c.upper[!c.id]
}
} else {
if(!drop.nugget){
c.names = c(c.names,"tau2")
c.lower = c(c.lower,0)
c.upper = c(c.upper,Inf)
}
}
drop.rhos = ifelse((!non.stat & !fit.nugget) | (non.stat & !fit.nugget),TRUE,FALSE)
if(any(grepl("(log.rhos|rhos)",c.names))){
if(drop.rhos){
c.id = grepl("(log.rhos|rhos)",c.names)
c.names = c.names[!c.id]
c.lower = c.lower[!c.id]
c.upper = c.upper[!c.id]
}
}
drop.kappa = ifelse(is.null(pars$fix.kappa),FALSE,ifelse(is.logical(pars$fix.kappa),ifelse(pars$fix.kappa,TRUE,FALSE),TRUE))
if(any(grepl("kappa",c.names))){
if(drop.kappa){
c.id = grepl("kappa",c.names)
c.names = c.names[!c.id]
c.lower = c.lower[!c.id]
c.upper = c.upper[!c.id]
}
}
names(c.lower) = c.names
names(c.upper) = c.names
list(lower = c.lower, upper = c.upper)
}
.gs = function(zV, XM, ZM, VM, hM, Qel, W,
use = list(penalty = "cs", ps.order = 2, Ztr = "RVA",
inverse = ".invEigen", tol = .Machine$double.neg.eps,
Amethod = "eigen", Adiag = TRUE),
pars = list(addbeta = FALSE, is.log.rhos = TRUE,
fix.nugget = NULL, fix.kappa = NULL,
model = "matern", nugget.tol = 1.0e-15),
method = list(use.reml = FALSE, use.profile = TRUE),
control = list(is.zero = 1.0e-6, is.inf = 10000, size.each.grid = 30,
size.sample = 40),
print.pars = FALSE){
if(is.null(use)) stop("List use must be specified!")
if(is.null(pars)) stop("List pars must be specified!")
if(is.null(method)) stop("List method must be specified!")
if(is.null(control)) stop("List control must be specified!")
if(is.null(use$penalty)) stop("Element penalty must be defined in use list!")
if(is.null(use$ps.order) & use$penalt=="ps") stop("Element ps.order must be defined in use list!")
if(is.null(use$Ztr)) stop("Element Ztr must be defined in use list!")
if(is.null(use$inverse)) stop("Element inverse must be defined in use list!")
if(is.null(use$tol)) stop("Element tol must be defined in use list!")
if(is.null(use$Amethod)) stop("Element Amethod must be defined in use list!")
if(is.null(use$Adiag)) stop("Element Adiag must be defined in use list!")
if(is.null(pars$addbeta)) stop("Element addbeta must be defined in pars list!")
if(is.null(pars$is.log.rhos)) stop("Element is.log.rhos must be defined in pars list!")
if(is.null(pars$model)) stop("Element model must be defined in pars list!")
if(is.null(pars$nugget.tol)) stop("Element nugget.tol must be defined in pars list!")
if(is.null(method$use.reml)) stop("Element use.reml must be defined in method list!")
if(is.null(method$use.profile)) stop("Element use.profile must be defined in method list!")
if(is.null(control$is.zero)) stop("Element is.zero must be defined in control list!")
if(is.null(control$is.inf)) stop("Element is.inf must be defined in control list!")
if(is.null(control$size.each.grid)) stop("Element size.each.grid must be defined in control list!")
if(is.null(control$size.sample)) stop("Element size.sample must be defined in control list!")
use$Adiag = ifelse(use$penalty=="ps", TRUE, use$Adiag)
c.theta = .getBounds(pars = pars, smooth.type=use$penalty)
c.npars = length(c.theta$lower)
c.nat = control$size.each.grid**c.npars
if(c.nat > 35**5){
c.base.size = 10
while(c.base.size**c.npars < 35**5){
c.base.size = c.base.size + 1
}
control$size.each.grid = c.base.size - 1
message(paste("The value size.each.grid**length(pars) cannot be greated than 50**5.\nChanging value to size.each.grid = ",c.base.size-1,".\nTo change manually modify control and pars lists!",sep=""))
}
c.start = c.theta
if(pars$addbeta){
c.start$lower = c(beta=rep(-Inf,ncol(XM)),c.theta$lower)
c.start$upper = c(beta=rep(+Inf,ncol(XM)),c.theta$upper)
}
c.init.pars = .randomParsGrid(c.theta, control = control)
c.ilL = NULL
for(i in 1:nrow(c.init.pars)){
message(paste("Step",i,"of",nrow(c.init.pars)))
aux = .tryCatchWE(.scpLoglik(
par = c.init.pars[i,],
beta.par = FALSE,
is.log.rhos = pars$is.log.rhos,
obs.z = zV,
dist.m = hM,
W = W,
Xs = XM,
Z = ZM,
V = VM,
fix.nugget = pars$fix.nugget,
fix.kappa = pars$fix.kappa,
pty = use$penalty,
Q.base = Qel,
Ztransf = use$Ztr,
model = pars$model,
ml = !method$use.reml,
use.profile = method$use.profile,
inverse.method = use$inverse,
itol = use$tol,
createA = use$Amethod,
diagA = use$Adiag,
nugget.tol = pars$nugget.tol,
print.pars = print.pars
))
c.ilL[i] = ifelse(!aux$error, aux$value, NA)
}
c.start$pars = c.init.pars[match(min(c.ilL, na.rm = TRUE),c.ilL),]
if(length(c.start$pars)==1)
names(c.start$pars) <- colnames(c.init.pars)
if(pars$addbeta){
c.beta = .betaFromPars(pars = c.start$pars, zV = zV, XM = XM, ZM = ZM, VM = VM, hM = hM, Qel = Qel, W,
use = use, pars.str = pars, use.reml = method$use.reml, only.beta = TRUE)
c.start$pars = c(c.beta,c.start$pars)
}
return(c.start)
}
.rightDistance = function(z,psi){
aux = cbind((z-psi)*(z>psi),-(z>psi))
colnames(aux) = paste(c("U(","V("),psi,")",sep="")
return(aux)
}
.matrixUV = function(z=NULL,psi){
zMatrix = matrix(0,nrow=length(z),ncol=2*length(psi))
cnamesz = rep("",ncol=2*length(psi))
for(i in 1:length(psi)){
aux = .rightDistance(z,psi[i])
cnamesz[(2*i-1):(2*i)] = colnames(aux)
zMatrix[,(2*i-1):(2*i)] = aux
}
colnames(zMatrix) = cnamesz
return(zMatrix)
}
.psiUpdate = function(beta, psi0, UVcols){
if(length(psi0)*2!=length(UVcols)) stop("Length of psi must match the number of columns of UV matrix!")
return(psi0 + (beta[seq(min(UVcols)+1,max(UVcols),2)])/(beta[seq(min(UVcols),max(UVcols)-1,2)]))
}
.psiIterative = function(cp.call, psi0, UVcols, obj, control.cp = list(maxiter = 20, maxskip = 10)){
psi0.v = NULL
UVcols.v = NULL
for(i in 1:length(psi0)){
psi0.v = c(psi0.v,psi0[[i]])
UVcols.v = c(UVcols.v,UVcols[[i]])
}
beta.step = matrix(0,nrow=control.cp$maxiter,ncol=ncol(obj$XM))
psi.step = matrix(0,nrow=control.cp$maxiter,ncol=length(psi0.v))
beta.step[1,] = rep(0,ncol(obj$XM))
psi.step[1,] = psi0.v
num.skip = 0
ii = 2
while(ii < (control.cp$maxiter + 1) && num.skip < control.cp$maxskip){
if(ii>2){
acc.psi = 0
newUV = NULL
for(i in 1:length(cp.call)){
mp = match("psi",names(cp.call[[i]]))
cp.call[[i]][[mp]] = psi.step[ii-1,(acc.psi+1):(acc.psi+length(psi0[[i]]))]
newUV[[i]] = eval(cp.call[[i]], parent.frame())
obj$XM[,UVcols[[i]]] = newUV[[i]]$UV
acc.psi = acc.psi + length(psi0[[i]])
}
}
cat(".")
beta.aux = .betaFromRSPZ(obj)
beta.step[ii,] = beta.aux$beta
psi.step[ii,] = .psiUpdate(beta.step[ii,], psi.step[ii-1,], UVcols.v)
ii = ii + 1
}
psi = psi.step[nrow(psi.step),]
acc.psi = 0
psi.l = NULL
for(i in 1:length(cp.call)){
psi.l[[i]] = psi[(acc.psi+1):(acc.psi+length(psi0[[i]]))]
acc.psi = acc.psi + length(psi0[[i]])
}
list(psi = psi.l, beta = beta.step[nrow(beta.step),], XM = obj$XM, UVcols = UVcols,
elements = beta.aux, history = list(beta = beta.step, psi = psi.step))
}
.RSPZfromPars = function(pars, zV, XM, ZM, VM, hM, Qel, W,
use = list(penalty = "cs", ps.order = 2, Ztr = "RVA",
inverse = ".invEigen", tol = .Machine$double.neg.eps,
Amethod = "eigen", Adiag = TRUE),
pars.str = list(fix.nugget = NULL, fix.kappa = NULL, is.log.rhos = TRUE,
model = "matern", nugget.tol = 1.0e-15)){
if(missing(Qel)) Qel <- NULL
if(missing(use)) use <- NULL
if(is.null(use)){
use = list(penalty = "none", inverse = ".invEigen", tol = .Machine$double.neg.eps)
}
if(is.null(use$penalty)) use$penalty <- "none"
if(is.null(use$intercept)) use$intercept <- TRUE
if(is.null(pars.str)){
stop("pars.str must be defined!")
} else {
bounds = .getBounds(pars = list(fix.nugget = pars.str$fix.nugget,
fix.kappa = pars.str$fix.kappa,
is.log.rhos = pars.str$is.log.rhos,
model = pars.str$model), smooth.type = use$penalty)
c.req = names(bounds$lower)
}
c.id.el = match(c.req,names(pars))
if(all(!is.na(c.id.el))){
pars = pars[c.id.el]
} else {
stop(paste("One or more required elements of pars were not specified!\nEnter ",paste(c.req[is.na(c.id.el)],sep="",collapse=", "),".",sep=""))
}
ac <- paste(formals(.autoCor)$model)[-1]
sv <- paste(formals(.semiVar)$model)[-1]
ul <- unique(c(ac,sv))
fl <- sv[is.na(match(sv,ac))]
non.stat <- any(grepl(paste("^",pars.str$model,"$",sep="",collapse=""),fl))
nugget.par <- ifelse(non.stat & is.null(pars.str$fix.nugget),TRUE,FALSE)
fit.nugget = ifelse(is.null(pars.str$fix.nugget),TRUE,FALSE)
fit.rhos = ifelse(fit.nugget,TRUE,FALSE)
fit.kappa = ifelse(is.null(pars.str$fix.kappa),TRUE,FALSE)
drop.nugget = ifelse(!non.stat | (non.stat & !fit.nugget),TRUE,FALSE)
drop.rhos = ifelse((!non.stat & !fit.nugget) | (non.stat & !fit.nugget),TRUE,FALSE)
if(use$penalty=="none") smooth <- FALSE else smooth <- TRUE
if(smooth){
if(use$penalty=="cs" | use$penalty=="ps"){
c.id.a = match(c("alpha1","alpha2"),names(pars))
} else if(use$penalty=="tps"){
c.id.a = match(c("alpha"),names(pars))
} else if(use$penalty!="none"){ stop("Non valid smooth type definition!") }
if(all(!is.na(c.id.a))) c.alpha.xy = pars[c.id.a] else stop("One or more elements of alpha not specified!")
c.alpha <- sum(c.alpha.xy)
c.rho.xy = c.alpha.xy/c.alpha
} else {
c.alpha.xy = NULL
c.alpha <- NULL
c.rho.xy = NULL
}
c.id.p = match("phi",names(pars))
if(all(!is.na(c.id.p))) c.phi = pars[c.id.p] else c.phi = NULL
if(is.null(pars.str$fix.kappa)){
c.id.k = grepl("kappa",names(pars))
if(any(c.id.k)) c.kappa = pars[c.id.k] else c.kappa = NULL
} else {
c.kappa = pars.str$fix.kappa
}
if(fit.nugget & !drop.nugget){
c.id.t = match("tau2",names(pars))
if(all(!is.na(c.id.t))) c.tau2 = pars[c.id.t] else c.tau2 = NULL
} else if(!is.null(pars.str$fix.nugget)){
c.tau2 <- NULL
c.rhos <- 1 - pars.str$fix.nugget
c.log.rhos <- log(c.rhos)
} else {
c.tau2 <- NULL
}
if(fit.rhos & !drop.rhos){
c.id.s = match(ifelse(pars.str$is.log.rhos,"log.rhos","rhos"),names(pars))
if(all(!is.na(c.id.s))){
if(pars.str$is.log.rhos){
c.log.rhos = pars[c.id.s]
c.rhos = exp(c.log.rhos)
} else {
c.rhos = pars[c.id.s]
c.log.rhos = log(c.rhos)
}
}
} else if(!is.null(pars.str$fix.nugget)){
c.tau2 <- NULL
c.rhos <- 1 - pars.str$fix.nugget
c.log.rhos <- log(c.rhos)
} else {
c.tau2 <- 0
c.rhos <- 1
c.log.rhos <- 0
}
if(!non.stat){
c.R.s = .createCov(dist.m = hM, phi = c.phi, sill = NULL, kappa = c.kappa, nugget = NULL, rhos = c.rhos,
model = pars.str$model, tol.nugget = pars.str$nugget.tol)
c.R.u = W%*%c.R.s%*%t(W)
} else {
c.sigma2s <- c.tau2/c.rhos
c.sigma2 <- c.sigma2s - c.tau2
c.Sigma.s = .createCov(dist.m = hM, phi = c.phi, sill = NULL, kappa = c.kappa, nugget = c.tau2, rhos = c.rhos,
model = pars.str$model, tol.nugget = pars.str$nugget.tol)
c.R.s = (1/c.sigma2s)%x%c.Sigma.s
c.Sigma.u = W%*%c.Sigma.s%*%t(W)
c.R.u = (1/c.sigma2s)%x%c.Sigma.u
}
c.H.u = .whichInverse(use$inverse, c.R.u, tol = use$tol, Re.values = TRUE)
if(smooth){
if(use$penalty=="tps"){
c.Ps <- Qel$K
c.l.Ps <- Qel
c.Z <- Qel$K
c.Gm = (1/c.alpha)%x%c.l.Ps$A
c.Gi = c.alpha%x%c.l.Ps$Ai
} else if(use$penalty=="cs" | use$penalty=="ps"){
c.Ps = .tensorPenalty(obj = Qel, alpha = c.rho.xy, type = "P")
c.l.Ps <- list(V = VM)
if(use$Ztr=="V" | use$Ztr=="VA"){
c.Z <- ZM
} else {
c.Z <- .generateZs(c.R.s, VM, itol = use$tol)
}
c.l.Ps$Ai <- .generateAi(c.Z, c.Ps, itol = use$tol)
c.l.Ps$A <- solve(c.l.Ps$Ai , tol = use$tol)
if(use$Ztr=="VA" | use$Ztr=="RVA"){
el.A <- eigen(c.l.Ps$A)
c.l.Ps$Asr <- el.A$vectors%*%diag(sqrt(el.A$values))%*%t(el.A$vectors)
c.Z <- c.Z%*%c.l.Ps$Asr
c.Gm = (1/c.alpha)%x%diag(1,nrow(c.l.Ps$A))
c.Gi = c.alpha%x%diag(1,nrow(c.l.Ps$A))
} else {
c.Gm = (1/c.alpha)%x%c.l.Ps$A
c.Gi = c.alpha%x%c.l.Ps$Ai
}
}
c.Zs <- W%*%c.Z
c.LL = c.Zs%*%c.Gm%*%t(c.Zs) + c.R.u
c.ZHZa = t(c.Zs)%*%c.H.u%*%c.Zs + c.Gi
c.ZHZai = .whichInverse(use$inverse, c.ZHZa, tol = use$tol, Re.values = TRUE)
c.tZ = c.ZHZai%*%t(c.Zs)%*%c.H.u
} else {
c.LL = c.R.u
c.ZHZa = NULL
c.ZHZai = NULL
c.tZ = NULL
}
c.In = diag(1,nrow(XM))
c.LLi = .whichInverse(use$inverse, c.LL, tol = use$tol, Re.values = TRUE)
list(alpha = if(exists("c.alpha.xy")){c.alpha.xy}, sigma2s=if(exists("c.sigma2s")){c.sigma2s}, sigma2s=if(exists("c.sigma2")){c.sigma2}, rhos = c.rhos, log.rhos = c.log.rhos, phi = c.phi, kappa = c.kappa, tau2 = c.tau2, R.s = c.R.s, R.u = c.R.u, H.u = c.H.u, Sigma.s = if(exists("c.Sigma.s")){c.Sigma.s}, Sigma.u = if(exists("c.Sigma.u")){c.Sigma.u}, W = W, XM = XM, Z = if(exists("c.Z")){c.Z}, Zs = if(exists("c.Zs")){c.Zs}, G = if(exists("c.Gm")){c.Gm}, Gi = if(exists("c.Gi")){c.Gi}, L = c.LL, Li = c.LLi, TZ = if(exists("c.tZ")){c.tZ}, ZHZa = if(exists("c.ZHZa")){c.ZHZa}, ZHZai = if(exists("c.ZHZai")){c.ZHZai}, zV = zV, use = use, pars.str = pars.str, penalty = if(use$penalty!="none"){list(P = c.Ps, decompose = c.l.Ps)}, fit.nugget = fit.nugget, Qel = Qel)
}
.betaFromPars = function(pars, zV, XM, ZM, VM, hM, Qel, W,
use = list(penalty = "cs", ps.order = 2, Ztr = "RVA",
inverse = ".invEigen", tol = .Machine$double.neg.eps,
Amethod = "eigen", Adiag = TRUE),
pars.str = list(fix.nugget = NULL, fix.kappa = NULL, is.log.rhos = TRUE,
model = "matern", nugget.tol = 1.0e-15),
use.reml = FALSE, only.beta = TRUE){
if(missing(Qel)) Qel <- NULL
if(missing(use)) use <- NULL
if(is.null(use)){
use = list(penalty = "none", inverse = ".invEigen", tol = .Machine$double.neg.eps)
}
if(is.null(use$penalty)) use$penalty <- "none"
if(is.null(use$intercept)) use$intercept <- TRUE
if(is.null(pars.str)){
stop("pars.str must be defined!")
} else {
bounds = .getBounds(pars = list(fix.nugget = pars.str$fix.nugget,
fix.kappa = pars.str$fix.kappa,
is.log.rhos = pars.str$is.log.rhos,
model = pars.str$model), smooth.type = use$penalty)
c.req = names(bounds$lower)
}
c.id.el = match(c.req,names(pars))
if(all(!is.na(c.id.el))){
pars = pars[c.id.el]
} else {
stop(paste("One or more required elements of pars were not specified!\nEnter ",paste(c.req[is.na(c.id.el)],sep="",collapse=", "),".",sep=""))
}
ac <- paste(formals(.autoCor)$model)[-1]
sv <- paste(formals(.semiVar)$model)[-1]
ul <- unique(c(ac,sv))
fl <- sv[is.na(match(sv,ac))]
non.stat <- any(grepl(paste("^",pars.str$model,"$",sep="",collapse=""),fl))
nugget.par <- ifelse(non.stat & is.null(pars.str$fix.nugget),TRUE,FALSE)
fit.nugget = ifelse(is.null(pars.str$fix.nugget),TRUE,FALSE)
fit.rhos = ifelse(fit.nugget,TRUE,FALSE)
fit.kappa = ifelse(is.null(pars.str$fix.kappa),TRUE,FALSE)
drop.nugget = ifelse(!non.stat | (non.stat & !fit.nugget),TRUE,FALSE)
drop.rhos = ifelse((!non.stat & !fit.nugget) | (non.stat & !fit.nugget),TRUE,FALSE)
if(use$penalty=="none") smooth <- FALSE else smooth <- TRUE
if(smooth){
if(use$penalty=="cs" | use$penalty=="ps"){
c.id.a = match(c("alpha1","alpha2"),names(pars))
} else if(use$penalty=="tps"){
c.id.a = match(c("alpha"),names(pars))
} else if(use$penalty!="none"){ stop("Non valid smooth type definition!") }
if(all(!is.na(c.id.a))) c.alpha.xy = pars[c.id.a] else stop("One or more elements of alpha not specified!")
c.alpha <- sum(c.alpha.xy)
c.rho.xy = c.alpha.xy/c.alpha
} else {
c.alpha.xy = NULL
c.alpha <- NULL
c.rho.xy = NULL
}
c.id.p = match("phi",names(pars))
if(all(!is.na(c.id.p))) c.phi = pars[c.id.p] else c.phi = NULL
c.fit.nugget <- fit.nugget
c.fit.kappa <- fit.kappa
if(is.null(pars.str$fix.kappa)){
c.id.k = grepl("kappa",names(pars))
if(any(c.id.k)) c.kappa = pars[c.id.k] else c.kappa = NULL
} else {
c.kappa = pars.str$fix.kappa
}
if(fit.nugget & !drop.nugget){
c.id.t = match("tau2",names(pars))
if(all(!is.na(c.id.t))) c.tau2 = pars[c.id.t] else c.tau2 = NULL
} else if(!is.null(pars.str$fix.nugget)){
c.tau2 <- NULL
c.rhos <- 1 - pars.str$fix.nugget
c.log.rhos <- log(c.rhos)
} else {
c.tau2 <- 0
c.rhos <- 1
c.log.rhos <- 0
}
if(fit.rhos & (!drop.rhos)){
c.id.s = match(ifelse(pars.str$is.log.rhos,"log.rhos","rhos"),names(pars))
if(all(!is.na(c.id.s))){
if(pars.str$is.log.rhos){
c.log.rhos = pars[c.id.s]
c.rhos = exp(c.log.rhos)
} else {
c.rhos = pars[c.id.s]
c.log.rhos = log(c.rhos)
}
}
} else if(!is.null(pars.str$fix.nugget)){
c.tau2 <- NULL
c.rhos <- 1 - pars.str$fix.nugget
c.log.rhos <- log(c.rhos)
} else {
c.tau2 <- NULL
c.rhos <- 1
c.log.rhos <- 0
}
if(!non.stat){
c.R.s = .createCov(dist.m = hM, phi = c.phi, sill = NULL, kappa = c.kappa, nugget = NULL, rhos = c.rhos,
model = pars.str$model, tol.nugget = pars.str$nugget.tol)
c.R.u = W%*%c.R.s%*%t(W)
} else {
c.sigma2s <- c.tau2/c.rhos
c.sigma2 <- c.sigma2s - c.tau2
c.Sigma.s = .createCov(dist.m = hM, phi = c.phi, sill = NULL, kappa = c.kappa, nugget = c.tau2, rhos = c.rhos,
model = pars.str$model, tol.nugget = pars.str$nugget.tol)
c.R.s = (1/c.sigma2s)%x%c.Sigma.s
c.Sigma.u = W%*%c.Sigma.s%*%t(W)
c.R.u = (1/c.sigma2s)%x%c.Sigma.u
}
c.H.u = .whichInverse(use$inverse, c.R.u, tol = use$tol, Re.values = TRUE)
c.In = diag(1,nrow(XM))
if(smooth){
if(use$penalty=="tps"){
c.Ps <- Qel$K
c.l.Ps <- Qel
c.Z <- Qel$K
c.Gm = (1/c.alpha)%x%c.l.Ps$A
c.Gi = c.alpha%x%c.l.Ps$Ai
} else if(use$penalty=="cs" | use$penalty=="ps"){
c.Ps = .tensorPenalty(obj = Qel, alpha = c.rho.xy, type = "P")
c.l.Ps <- list(V = VM)
if(use$Ztr=="V" | use$Ztr=="VA"){
c.Z <- ZM
} else {
c.Z <- .generateZs(c.R.s, VM, itol = use$tol)
}
c.l.Ps$Ai <- .generateAi(c.Z, c.Ps, itol = use$tol)
c.l.Ps$A <- solve(c.l.Ps$Ai , tol = use$tol)
if(use$Ztr=="VA" | use$Ztr=="RVA"){
el.A <- eigen(c.l.Ps$A)
c.l.Ps$Asr <- el.A$vectors%*%diag(sqrt(el.A$values))%*%t(el.A$vectors)
c.Z <- c.Z%*%c.l.Ps$Asr
c.Gm = (1/c.alpha)%x%diag(1,nrow(c.l.Ps$A))
c.Gi = c.alpha%x%diag(1,nrow(c.l.Ps$A))
} else {
c.Gm = (1/c.alpha)%x%c.l.Ps$A
c.Gi = c.alpha%x%c.l.Ps$Ai
}
}
c.Zs <- W%*%c.Z
c.LL = c.Zs%*%c.Gm%*%t(c.Zs) + c.R.u
c.ZHZa = t(c.Zs)%*%c.H.u%*%c.Zs + c.Gi
c.ZHZai = .whichInverse(use$inverse, c.ZHZa, tol = use$tol, Re.values = TRUE)
c.tZ = c.ZHZai%*%t(c.Zs)%*%c.H.u
c.pZ = c.Zs%*%c.tZ
c.IpZ = c.In - c.pZ
} else {
c.Zs <- NULL
c.LL = c.R.u
c.ZHZa = NULL
c.ZHZai = NULL
c.tZ = NULL
c.pZ = NULL
c.IpZ = c.In
}
c.LLi = .whichInverse(use$inverse, c.LL, tol = use$tol, Re.values = TRUE)
c.XLX = t(XM)%*%c.LLi%*%XM
c.XLXi = .whichInverse(use$inverse,c.XLX, tol = use$tol, Re.values = TRUE)
c.tXs = c.XLXi%*%t(XM)%*%c.LLi
c.beta.s = c.tXs%*%zV
names(c.beta.s) = paste("beta",1:length(c.beta.s),sep="")
if(only.beta){
names(c.beta.s) <- colnames(XM)
return(c.beta.s)
} else {
if(smooth){
c.r = c.tZ%*%(zV-XM%*%c.beta.s)
c.yf = XM%*%c.beta.s + c.Zs%*%c.r
} else {
c.r <- NULL
c.yf = XM%*%c.beta.s
}
if(use$is.X!="none"){
c.Vdim <- ifelse(use$is.X=="tensor",3+ifelse(use$intercept,1,0),2+ifelse(use$intercept,1,0))
if(ncol(XM)>c.Vdim){
c.C = XM[,-((ncol(XM)-(c.Vdim-1)):ncol(XM))]
c.X = XM[,(ncol(XM)-(c.Vdim-1)):ncol(XM)]
if(smooth){
c.pZ = c.Zs%*%c.tZ
c.IpZ = c.In - c.pZ
} else {
c.pZ = NULL
c.IpZ = c.In
}
c.XHX = t(c.X)%*%c.H.u%*%c.IpZ%*%c.X
c.XHXi = .whichInverse(use$inverse, c.XHX, tol = use$tol, Re.values = TRUE)
c.tX = c.XHXi%*%t(c.X)%*%c.H.u%*%c.IpZ
c.pX = c.X%*%c.tX
c.IpX <- c.In - c.pX
if(smooth){
c.Sxz = c.pX + c.pZ%*%c.IpX
} else {
c.Sxz <- c.pX
}
c.ISxz = c.In - c.Sxz
c.CHISC = t(c.C)%*%c.H.u%*%c.ISxz%*%c.C
c.CHISCi = .whichInverse(use$inverse, c.CHISC, tol = use$tol, Re.values = TRUE)
c.tC = c.CHISCi%*%t(c.C)%*%c.H.u%*%c.ISxz
c.pC = c.C%*%c.tC
c.IpC = c.In-c.pC
c.Mxz = c.pC + c.Sxz%*%c.IpC
if(smooth){
c.gf = c.X%*%(c.beta.s[(ncol(XM)-(c.Vdim-1)):ncol(XM)]) + c.Zs%*%c.r
} else {
c.gf = c.X%*%(c.beta.s[(ncol(XM)-(c.Vdim-1)):ncol(XM)])
}
} else {
c.X <- XM
if(smooth){
c.pZ = c.Zs%*%c.tZ
c.IpZ = c.In - c.pZ
} else {
c.pZ = NULL
c.IpZ = c.In
}
c.XHX = t(c.X)%*%c.H.u%*%c.IpZ%*%c.X
c.XHXi = .whichInverse(use$inverse, c.XHX, tol = use$tol, Re.values = TRUE)
c.tX = c.XHXi%*%t(c.X)%*%c.H.u%*%c.IpZ
c.pX = c.X%*%c.tX
c.IpX <- c.In - c.pX
if(smooth){
c.Sxz = c.pX + c.pZ%*%c.IpX
} else {
c.Sxz <- c.pX
}
c.Mxz = c.Sxz
c.ISxz = c.In - c.Sxz
c.gf = c.yf
}
} else {
c.Mxz = XM%*%c.tXs
c.Sxz = NULL
c.gf = NULL
}
c.epsilon = zV - XM%*%c.beta.s
c.RSS = t(c.epsilon)%*%c.LLi%*%(c.epsilon)
if(!use.reml) c.df = nrow(XM) else c.df = (nrow(XM)-ncol(XM))
c.sigma2s = c.RSS/c.df
c.tau2 <- (1-c.rhos)*c.sigma2s
c.sigma2 = c.sigma2s - c.tau2
c.Sigma.u = c.sigma2s%x%c.R.u
colnames(c.XLX) <- colnames(XM)
rownames(c.XLX) <- colnames(XM)
colnames(c.XLXi) <- colnames(XM)
rownames(c.XLXi) <- colnames(XM)
names(c.beta.s) <- colnames(XM)
list(alpha = if(exists("c.alpha.xy")){c.alpha.xy}, sigma2s = c.sigma2s, sigma2 = c.sigma2, rhos = c.rhos, log.rhos = c.log.rhos, phi = c.phi, kappa = c.kappa, tau2 = c.tau2, beta = c.beta.s, r = c.r, R = c.R.u, g = if(exists("c.gf")){c.gf}, fitted.values = if(exists("c.yf")){c.yf}, S = if(exists("c.Sxz")){c.Sxz}, M = if(exists("c.Mxz")){c.Mxz}, X = XM, Z = if(exists("c.Z")){c.Z}, Zs = if(exists("c.Zs")){c.Zs}, W = W, XLX = c.XLX, XLXi = c.XLXi, var.beta = c.sigma2s%x%c.XLXi, conf.error = if(!only.beta){sqrt(diag(c.Mxz%*%c.Sigma.u%*%t(c.Mxz)))}, pred.error = if(!only.beta){sqrt(diag(c.Sigma.u + c.Mxz%*%c.Sigma.u%*%t(c.Mxz)))}, penalty = if(use$penalty!="none"){list(P = c.Ps, decompose = c.l.Ps)})
}
}
.betaFromRSPZ = function(obj){
c.In = diag(1,nrow(obj$XM))
c.XLX = t(obj$XM)%*%obj$Li%*%obj$XM
c.XLXi = .whichInverse(obj$use$inverse, c.XLX, tol = obj$use$tol, Re.values = TRUE)
c.tXs = c.XLXi%*%t(obj$XM)%*%obj$Li
c.beta.s = c.tXs%*%obj$zV
names(c.beta.s) = paste("beta",1:length(c.beta.s),sep="")
list(beta = c.beta.s, XLX = c.XLX, XLXi = c.XLXi, TX = c.tXs)
}
.fitFromRSPZ = function(RZ, IP, use.reml = FALSE){
c.fit.nugget = RZ$fit.nugget
c.Ps = RZ$penalty$P
c.l.Ps = RZ$penalty$decompose
use = RZ$use
if(use$penalty=="none") smooth <- FALSE else smooth <- TRUE
pars.str = RZ$pars.str
c.rhos = c(rhos = RZ$rhos)
c.log.rhos = c(log.rhos = RZ$log.rhos)
c.tau2 = RZ$tau2
XM = IP$XM
zV = RZ$zV
c.LLi = RZ$Li
c.XLX = t(XM)%*%c.LLi%*%XM
c.XLXi = .whichInverse(use$inverse, c.XLX, tol = use$tol, Re.values = TRUE)
c.R.u = RZ$R.u
c.H.u = RZ$H.u
c.Sigma.u = RZ$Sigma.u
c.Z = RZ$Z
c.Zs = RZ$Zs
c.tZ = RZ$TZ
c.tXs = c.XLXi%*%t(XM)%*%c.LLi
c.beta.s = IP$beta
c.In = diag(1,nrow(XM))
if(smooth){
c.r = c.tZ%*%(zV-XM%*%c.beta.s)
c.yf = XM%*%c.beta.s + c.Zs%*%c.r
} else {
c.r <- NULL
c.yf = XM%*%c.beta.s
}
if(use$is.X!="none"){
c.Vdim <- ifelse(use$is.X=="tensor",3+ifelse(use$intercept,1,0),2+ifelse(use$intercept,1,0))
if(ncol(XM)>c.Vdim){
c.C = XM[,-((ncol(XM)-(c.Vdim-1)):ncol(XM))]
c.X = XM[,(ncol(XM)-(c.Vdim-1)):ncol(XM)]
if(smooth){
c.pZ = c.Zs%*%c.tZ
c.IpZ = c.In - c.pZ
} else {
c.pZ = NULL
c.IpZ = c.In
}
c.XHX = t(c.X)%*%c.H.u%*%c.IpZ%*%c.X
c.XHXi = .whichInverse(use$inverse, c.XHX, tol = use$tol, Re.values = TRUE)
c.tX = c.XHXi%*%t(c.X)%*%c.H.u%*%c.IpZ
c.pX = c.X%*%c.tX
c.IpX <- c.In - c.pX
if(smooth){
c.Sxz = c.pX + c.pZ%*%c.IpX
} else {
c.Sxz <- c.pX
}
c.ISxz = c.In - c.Sxz
c.CHISC = t(c.C)%*%c.H.u%*%c.ISxz%*%c.C
c.CHISCi = .whichInverse(use$inverse, c.CHISC, tol = use$tol, Re.values = TRUE)
c.tC = c.CHISCi%*%t(c.C)%*%c.H.u%*%c.ISxz
c.pC = c.C%*%c.tC
c.IpC = c.In-c.pC
c.Mxz = c.pC + c.Sxz%*%c.IpC
if(smooth){
c.gf = c.X%*%(c.beta.s[(ncol(XM)-(c.Vdim-1)):ncol(XM)]) + c.Zs%*%c.r
} else {
c.gf = c.X%*%(c.beta.s[(ncol(XM)-(c.Vdim-1)):ncol(XM)])
}
} else {
c.X <- XM
if(smooth){
c.pZ = c.Zs%*%c.tZ
c.IpZ = c.In - c.pZ
} else {
c.pZ = NULL
c.IpZ = c.In
}
c.XHX = t(c.X)%*%c.H.u%*%c.IpZ%*%c.X
c.XHXi = .whichInverse(use$inverse, c.XHX, tol = use$tol, Re.values = TRUE)
c.tX = c.XHXi%*%t(c.X)%*%c.H.u%*%c.IpZ
c.pX = c.X%*%c.tX
c.IpX <- c.In - c.pX
if(smooth){
c.Sxz = c.pX + c.pZ%*%c.IpX
} else {
c.Sxz <- c.pX
}
c.Mxz = c.Sxz
c.ISxz = c.In - c.Sxz
c.gf = c.yf
}
} else {
c.Mxz = XM%*%c.tXs
c.Sxz = NULL
c.gf = NULL
}
c.epsilon = zV - XM%*%c.beta.s
c.RSS = t(c.epsilon)%*%c.LLi%*%(c.epsilon)
if(!use.reml) c.df = nrow(XM) else c.df = (nrow(XM)-ncol(XM))
c.sigma2s = c.RSS/c.df
c.tau2 <- (1-c.rhos)*c.sigma2s
c.sigma2 = c.sigma2s - c.tau2
c.Sigma.u = c.sigma2s%x%c.R.u
colnames(c.XLX) <- colnames(XM)
rownames(c.XLX) <- colnames(XM)
colnames(c.XLXi) <- colnames(XM)
rownames(c.XLXi) <- colnames(XM)
names(c.beta.s) <- colnames(XM)
list(alpha = RZ$alpha, sigma2s = c.sigma2s, sigma2 = c.sigma2, phi = RZ$phi, kappa = RZ$kappa, tau2 = c.tau2, rhos = c.rhos, log.rhos = c.log.rhos, beta = c.beta.s, r = c.r, psi = IP$psi, history.psi = IP$history, R = c.R.u, g = if(exists("c.gf")){c.gf}, fitted.values = if(exists("c.yf")){c.yf}, S = if(exists("c.Sxz")){c.Sxz}, M = if(exists("c.Mxz")){c.Mxz}, X = XM, Z = if(exists("c.Z")){c.Z}, Zs = if(exists("c.Zs")){c.Zs}, W = RZ$W, XLX = c.XLX, XLXi = c.XLXi, var.beta = c.sigma2s%x%c.XLXi, conf.error = sqrt(diag(c.Mxz%*%c.Sigma.u%*%t(c.Mxz))), pred.error = sqrt(diag(c.Sigma.u + c.Mxz%*%c.Sigma.u%*%t(c.Mxz))), penalty = if(use$penalty!="none"){list(P = c.Ps, decompose = c.l.Ps)})
}
.scpInitialize = function(zV, XM, ZM, VM, hM, Qel, W, ch, use, pars, initial, method,
control.optim, control.gs, control.ch, print.pars){
if(missing(use)){
use = list(penalty = "none", is.X = "none", ps.order = 2, Ztr = "RV",
inverse = "solve", tol = .Machine$double.neg.eps,
Amethod = "eigen", Adiag = TRUE)
} else {
if(is.null(use)){
stop("List use must be specified!")
} else {
if(is.null(use$penalty)){
use$penaly = "none"
}
if(use$penalty=="cs" | use$penalty=="ps"){
use$is.X = "tensor"
}
if(use$penalty=="tps"){
use$is.X = "tps"
}
if(is.null(use$is.X) & use$penalty=="none"){
use$is.X = "none"
}
if(is.null(use$ps.order) & use$penalty=="ps"){
use$ps.order = 2
}
if(is.null(use$Ztr)){
use$Ztr = "RV"
}
if(is.null(use$inverse)){
use$inverse = "solve"
}
if(is.null(use$tol)){
use$tol = .Machine$double.neg.eps
}
if(is.null(use$Amethod)){
use$Amethod = "eigen"
}
if(is.null(use$Adiag)){
use$Adiag = ifelse(use$penalty=="ps", TRUE, TRUE)
}
}
}
if(missing(pars)){
pars = list(fix.nugget = NULL, fix.kappa = NULL, is.log.rhos = TRUE,
model = "matern", nugget.tol = 1.0e-15)
} else {
if(is.null(pars)){
stop("List pars must be specified!")
} else {
if(is.null(pars$model)){
pars$model = "matern"
}
if(is.null(pars$nugget.tol)){
pars$nugget.tol = 1.0e-15
}
}
}
use.gs = FALSE
alt.gs = FALSE
if(missing(initial)){
use.gs = TRUE
initial = NULL
}
if(is.null(initial)) message("Initial values not specified. Using internal search!")
if(missing(method)){
method = list(use.reml = FALSE, use.profile = TRUE, linearize = "eigen")
} else {
if(is.null(method)){
stop("List method must be specified!")
} else {
if(is.null(method$use.reml)){
method$use.reml = FALSE
}
if(is.null(method$use.profile)){
method$use.profile = TRUE
}
if(is.null(method$linearize)){
method$linearize = "eigen"
}
}
}
if(missing(control.optim)){
control.optim = list(trace = 5, maxit = 300)
} else {
if(is.null(control.optim)){
control.optim = list()
}
}
if(missing(control.gs)){
control.gs = list(is.zero = 1.0e-6, is.inf = 10000, size.each.grid = 30, size.sample = 40, num.attempts = 5)
} else {
if(is.null(control.gs)){
stop("List control.gs must be specified!")
} else {
if(is.null(control.gs$is.zero)){
control.gs$is.zero = 1.0e-6
}
if(is.null(control.gs$is.inf)){
control.gs$is.inf = 10000
}
if(is.null(control.gs$size.each.grid)){
control.gs$size.each.grid = 30
}
if(is.null(control.gs$size.sample)){
control.gs$size.sample = 40
}
if(is.null(control.gs$num.attempts)){
control.gs$num.attempts = 5
}
if(is.null(control.gs$replace.optim)){
control.gs$replace.optim = FALSE
}
}
}
if(missing(ch)){
control.ch = list(maxiter = 20, maxskip = 10)
} else {
if(is.null(control.ch)){
stop("List control.ch must be specified!")
} else {
if(is.null(control.ch$maxiter)){
control.ch$maxiter = 20
}
if(is.null(control.ch$maxskip)){
control.ch$maxskip = 10
}
}
}
if(missing(print.pars)){ print.pars = FALSE }
use$Adiag = ifelse(use$penalty=="ps", TRUE, use$Adiag)
cycle <- 1
while(cycle <= control.gs$num.attempts){
cat(".")
message(paste("Computing: cycle",cycle))
c.initial = unlist(initial)
c.theta = .getBounds(pars = pars, smooth.type = use$penalty)
c.req = names(c.theta$lower)
c.start = c.theta
c.id.el = match(c.req,names(c.initial))
if(all(!is.na(c.id.el))){
c.start$pars = c.initial[c.id.el]
} else {
use.gs = TRUE
c.start$pars = c.initial[!is.na(match(names(c.initial),c.req))]
}
if(!is.null(initial)){
if(any(grepl("beta",names(initial)))){
if(all(!is.na(match(names(initial$beta),paste("beta",1:ncol(XM),sep=""))))){
pars$addbeta = TRUE
c.start$pars = c(initial$beta,c.start$pars)
c.start$lower = c(beta=rep(-Inf,ncol(XM)),c.start$lower)
c.start$upper = c(beta=rep(+Inf,ncol(XM)),c.start$upper)
} else {
message("Length of starting values for beta does not match the required!\nUsing internal search for starting values not given!")
pars$addbeta = TRUE
use.gs = TRUE
}
} else {
pars$addbeta = FALSE
}
} else {
pars$addbeta = FALSE
}
if(use.gs | control.gs$replace.optim){
message("One or more starting values were not provided.\nUsing internal search!")
c.gstart = .gs(zV = zV, XM = XM, ZM = ZM, VM = VM, hM = hM, Qel = Qel, W = W,
use = use, pars = pars,
method = method, control = control.gs,
print.pars = FALSE)
c.gs = c.gstart
if(use.gs){
c.id.given = match(names(c.gstart$pars),names(c.start$pars))
if(sum(!is.na(c.id.given))>0){
c.gstart$pars[!is.na(c.id.given)] = c.start$pars[c.id.given[!is.na(c.id.given)]]
} else {
alt.gs = TRUE
}
c.start = c.gstart
}
}
message(paste("Estimating ",paste(names(c.start$pars),sep="",collapse=", "),sep="",collapse=""))
message(paste("Starting at ",paste(round(c.start$pars,5),sep="",collapse=", "),sep="",collapse=""))
requireNamespace("stats", quietly = TRUE)
message(paste(ifelse(method$use.reml,"REML","ML"),"estimation.",ifelse(method$use.profile,"Use profiling.\n","Not profiling.\n"),sep=" "))
aux = .tryCatchWE(
stats::optim(
par = c.start$pars,
.scpLoglik,
gr = NULL,
beta.par = pars$addbeta,
is.log.rhos = pars$is.log.rhos,
obs.z = zV,
dist.m = hM,
W = W,
Xs = XM,
Z = ZM,
V = VM,
fix.nugget = pars$fix.nugget,
fix.kappa = pars$fix.kappa,
pty = use$penalty,
Q.base = Qel,
Ztransf = use$Ztr,
model = pars$model,
ml = !method$use.reml,
use.profile = method$use.profile,
inverse.method = use$inverse,
itol = use$tol,
createA = use$Amethod,
diagA = use$Adiag,
nugget.tol = pars$nugget.tol,
print.pars = print.pars,
method = "L-BFGS-B", control = control.optim,
hessian = TRUE, lower = c.start$lower, upper = c.start$upper
)
)
cycle <- cycle + 1
if(!aux$error){
if(aux$value$convergence==0){
cycle <- control.gs$num.attempts + 1
}
}
}
if(!aux$error){
c.fit = aux$value
} else {
if((use.gs & alt.gs) | control.gs$replace.optim){
c.fit = list()
if(use.gs & alt.gs){
c.fit$par = c.start$pars
}
if(control.gs$replace.optim){
c.fit$par = c.gs$pars
}
message("Computing the Hessian")
c.fit$hessian = stats::optimHess(
par = c.fit$par,
.scpLoglik,
gr = NULL,
beta.par = pars$addbeta,
is.log.rhos = pars$is.log.rhos,
obs.z = zV,
dist.m = hM,
W = W,
Xs = XM,
Z = ZM,
V = VM,
fix.nugget = pars$fix.nugget,
fix.kappa = pars$fix.kappa,
pty = use$penalty,
Q.base = Qel,
Ztransf = use$Ztr,
model = pars$model,
ml = !method$use.reml,
use.profile = method$use.profile,
inverse.method = use$inverse,
itol = use$tol,
createA = use$Amethod,
diagA = use$Adiag,
nugget.tol = pars$nugget.tol,
print.pars = print.pars,
control = control.optim
)
} else {
message("Error in optim procedure!")
stop(aux$value)
}
}
if(is.null(ch)){
message("Obtaining mean estimates")
c.out = .betaFromPars(pars = c.fit$par, zV = zV, XM = XM, ZM = ZM, VM = VM, hM = hM, Qel = Qel, W = W,
use = use, pars.str = pars, use.reml = method$use.reml, only.beta = FALSE)
} else {
message("Obtaining mean estimates")
c.elem = .RSPZfromPars(pars = c.fit$par, zV = zV, XM = XM, ZM = ZM, VM = VM, hM = hM, Qel = Qel, W = W,
use = use, pars.str = pars)
message("Obtaining change-points")
c.psi.l = .psiIterative(cp.call = ch$cp.call, psi0 = ch$psi0, UVcols = ch$UVcols,
obj = c.elem, control.cp = control.ch)
message("Obtaining fitted model")
c.out = .fitFromRSPZ(RZ = c.elem, IP = c.psi.l, use.reml = method$use.reml)
}
if(!aux$error){
c.out$optim = c.fit
} else {
c.out$gs = c.fit
}
cat("\n")
return(c.out)
}
landim1 <- load(ifelse(file.exists("data/landim1.rda"),"data/landim1.rda",file.path(getwd(), "data/landim1.rda")))
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.