Skip to content

Commit

Permalink
mgks now based on distances and treated as "si" when inner == TRUE
Browse files Browse the repository at this point in the history
  • Loading branch information
Matteo Fasiolo committed Mar 29, 2024
1 parent 3dce73f commit 82f48b6
Show file tree
Hide file tree
Showing 2 changed files with 64 additions and 53 deletions.
115 changes: 63 additions & 52 deletions R/I_prepareInnerNested.R
Expand Up @@ -13,25 +13,26 @@

da <- length( alpha )
prange <- (sm$first.para:sm$last.para)[1:da]
Va <- gObj$Vp[prange, prange, drop = FALSE]

type <- class(o)[1]
if(type == "si"){
alpha <- drop(B %*% alpha)
Va <- B %*% Va %*% t(B)
se <- sqrt(pmax(0, diag(Va)))
edf <- sum(gObj$edf[prange])
ylabel <- .subEDF(paste0("proj_coef(", sm$term, ")"), edf)
xlabel <- "Index"
out <- list("fit" = alpha, "x" = 1:da, "se" = se,
xlab = xlabel, ylab = ylabel, main = NULL, type = "si")

# Remove scale parameter of inner transformation
if( type != "si" ){
prange <- prange[-1]
da <- da-1
alpha <- alpha[-1]
if(is.null(B)){
B <- diag(1, nrow = length(alpha))
}
}

Va <- gObj$Vp[prange, prange, drop = FALSE]

if( type == "nexpsm" ){
inner <- expsmooth(y = si$x, Xi = si$X, beta = alpha[-1], deriv = 1)
inner <- expsmooth(y = si$x, Xi = si$X, beta = alpha, deriv = 1)
fit <- inner$d0
Jac <- inner$d1
nobs <- length(fit)
se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
se <- sqrt(pmax(0, rowSums((Jac %*% Va) * Jac)))
edf <- sum(gObj$edf[prange[-1]])
ylabel <- .subEDF(paste0("expsm(", sm$term, ")"), edf)
xlabel <- "Index"
Expand All @@ -49,46 +50,56 @@
"p.resid" = si$x[ii], "raw" = ii,
"xlim" = xlim,
xlab = xlabel, ylab = ylabel, main = NULL, type = "nexpsm")
} else {
alpha <- drop(B %*% alpha)
Va <- B %*% Va %*% t(B)
se <- sqrt(pmax(0, diag(Va)))
edf <- sum(gObj$edf[prange])
ylabel <- .subEDF(paste0("Inner_coef(", sm$term, ")"), edf)
xlabel <- "Index"
out <- list("fit" = alpha, "x" = 1:da, "se" = se,
xlab = xlabel, ylab = ylabel, main = NULL, type = "si")
}
if( type == "mgks" ){
d <- ncol(si$X0)
if( d != 2 ){ # ONLY 2D case handled at the moment!!
return( NULL )
}
if( !is.null(xlim) ) {
xlim <- sort(xlim)
} else {
xlim <- range(si$X[ , 1])
}
if( !is.null(ylim) ) {
ylim <- sort(ylim)
} else {
ylim <- range(si$X[ , 2])
}

xx <- rep(seq(xlim[1], xlim[2], length.out = n), n)
yy <- rep(seq(ylim[1], ylim[2], length.out = n), rep(n, n))

X <- cbind(xx, yy)
si$x <- as.matrix(si$x)
if( ncol(si$x) > 1 ){
si$x <- colMeans(si$x)
}

inner <- mgks(y = si$x, X = X, X0 = si$X0, beta = alpha[-1], deriv = 1)
fit <- inner$d0
Jac <- inner$d1
se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
edf <- sum(gObj$edf[prange[-1]])

mainlab <- .subEDF(paste0("mgks(", sm$term, ")"), edf)
ylabel <- "X[ , 2]"
xlabel <- "X[ , 1]"
out <- list("fit" = fit, "X" = si$X, "se" = se, "x" = xx, "y" = yy,
"p.resid" = si$x, "X0" = si$X0,
"xlim" = xlim, "ylim" = ylim,
"xlab" = xlabel, "ylab" = ylabel, "main" = mainlab, type = "mgks")
}
# NOT CLEAR HOW TO DO THIS WITH DISTANCEs
# if( type == "mgks" ){
# d <- ncol(si$X0)
# if( d != 2 ){ # ONLY 2D case handled at the moment!!
# return( NULL )
# }
# if( !is.null(xlim) ) {
# xlim <- sort(xlim)
# } else {
# xlim <- range(si$X[ , 1])
# }
# if( !is.null(ylim) ) {
# ylim <- sort(ylim)
# } else {
# ylim <- range(si$X[ , 2])
# }
#
# xx <- rep(seq(xlim[1], xlim[2], length.out = n), n)
# yy <- rep(seq(ylim[1], ylim[2], length.out = n), rep(n, n))
#
# X <- cbind(xx, yy)
# si$x <- as.matrix(si$x)
# if( ncol(si$x) > 1 ){
# si$x <- colMeans(si$x)
# }
#
# inner <- mgks(y = si$x, X = X, X0 = si$X0, beta = alpha[-1], deriv = 1)
# fit <- inner$d0
# Jac <- inner$d1
# se <- sqrt(pmax(0, rowSums((Jac %*% Va[-1, -1, drop = FALSE]) * Jac)))
# edf <- sum(gObj$edf[prange[-1]])
#
# mainlab <- .subEDF(paste0("mgks(", sm$term, ")"), edf)
# ylabel <- "X[ , 2]"
# xlabel <- "X[ , 1]"
# out <- list("fit" = fit, "X" = si$X, "se" = se, "x" = xx, "y" = yy,
# "p.resid" = si$x, "X0" = si$X0,
# "xlim" = xlim, "ylim" = ylim,
# "xlab" = xlabel, "ylab" = ylabel, "main" = mainlab, type = "mgks")
# }

return(out)

Expand Down
2 changes: 1 addition & 1 deletion R/I_prepareOuterNested.R
Expand Up @@ -23,7 +23,7 @@
trnam <- "expsm"
}
if( type == "mgks" ){
raw <- mgks(y = si$x, X = si$X, X0 = si$X0, beta = alpha[-1])$d0
raw <- mgks(y = si$x, dist = si$dist, beta = alpha[-1])$d0
trnam <- "mgks"
}

Expand Down

0 comments on commit 82f48b6

Please sign in to comment.