R/KnePosSarEstimation_R1.R

In christophrust/FunRegPoI: Functional Linear Regression with Points of Impact

KnePosSarEstimation_R1 <-
function(Y, X_mat, add.vars, potPoI, searchMethod, grd, dom, plotting = FALSE,
         maxK = 40, maxPoI=10, efSelection = NULL, intercept = FALSE){
    ## #####################################################################
    ## Function performing the R1-estimate using the Eigenfunction-Expansion


    ## todo (additional variables currently only implemented in CraKneSaPoI):
    add.var.coef <- NULL
    
    ## Calculate centered Y and centered X
    p     <- nrow(X_mat)
    Y_c   <- scale(Y, scale = FALSE)
    X_c   <- scale(t(X_mat), scale = FALSE) # matplot(x = grd, y = X_c, type = "l"); dim(X_c); dim(X_mat)

    ## Calculate X_i(potPoI_j)
    PoIXValues <- scale(t(X_mat))[ , potPoI , drop = F]


    ## Choose from potPoI the PoIChoice searchMethod <- "dirSearch" searchMethod <- "fullSearch"
    PoIChoice <- if (searchMethod == "dirSearch"){
                     selectPoI_dirSearch(Y=scale(Y), X_mat=NULL, PoIXValues = PoIXValues, 
                                         potPoI = potPoI, k = 0, K = 0,
                                         maxPoI = maxPoI, nbest = 1, intercept = intercept, plotting = plotting)
                 } else {
                     selectPoI_fullSearch(Y=scale(Y), X_mat=NULL, PoIXValues = PoIXValues, 
                                          potPoI = potPoI, k = 0, K = 0,
                                          maxPoI = maxPoI, nbest = 1, intercept = intercept, plotting = plotting)
                 }
    S_choice  <- length(PoIChoice)


    ## Eigenfunctions, Basis Coefficients and Eigenvalues via PCA
    if (is.null(efSelection)){
        efSelection  <- selectEFwithCumPropVar(X_mat = X_mat)
    }
    scores       <- X_c %*% efSelection$evecs * 1/p

    ## Estimate Beta(t) and Beta_i for PoI_i
    ## select number of eigenfunctions wrt to gcv
    estDataFrame <- createEstDFforDirectSearch(Y = Y_c, X = scores[ ,1:maxK], PoI_vals = X_c[ ,PoIChoice, drop=FALSE])

    
    gcvs <- vapply(0:maxK , function(k){
        lmEst <- lm.fit(y = estDataFrame[,1], x = as.matrix(estDataFrame[,c(2:(k+1),(maxK+2):ncol(estDataFrame))]))
        calcGCVlm(lmEst)[2]
    },0.5)
    
    K_choice <- which.min(gcvs)-1

    
    ## get the final model estimates names(estDataFrame)
    if(K_choice != 0){
        estDataFrame <- createEstDFforDirectSearch(Y = Y_c, X = scores[ , 1:K_choice, drop=FALSE], PoI_vals = X_c[ , PoIChoice, drop=FALSE])
    } else {
        estDataFrame <- createEstDFforDirectSearch(Y = Y_c, X = NULL, PoI_vals = X_c[ , PoIChoice, drop=FALSE])
    }
    
    ## estimate Y againgst eigenfunctions and PoI choices
    LMandBIC    <- calcPoIandScoreLM(estDataFrame, s = S_choice, k = K_choice, K = K_choice, S = maxPoI)
    lm.fit      <- LMandBIC$lmFit
    coefnames   <- names(lm.fit$coeffic)
    score_names <- coefnames[substr(coefnames , 1,3)!="PoI"]
    poi_names   <- coefnames[substr(coefnames , 1,3)=="PoI"]
    estEF       <- efSelection$evecs[ , score_names]
    estPsi      <- lm.fit$coefficients[score_names]
    estPoI      <- lm.fit$coefficients[poi_names]
    
    ## Calculate beta_hat(t), PoI_hat, Y_hat
    beta_hat   <- calcScoreBetaFun(ef = estEF, coefs = estPsi, p=p)
    if (S_choice > 0) {
        ## Y_i = integral X_i(t) beta_hat(t) dt + sum_k beta_hat_k * X_i(tau_hat_k)
        Y_hat  <- (X_c %*% beta_hat(grd))/ p + X_c[, PoIChoice , drop=FALSE] %*% estPoI 
        
    } else {
        Y_hat  <- (X_c %*% beta_hat(grd))/ p
    }
    hat_matrix <- NULL
    list(selS = S_choice, 
         K_choice = K_choice,
         estPoI = estPoI,
         betaAddVar = add.var.coef,
         estTau = ind2grd(PoIChoice , grd),
         estTauGrd = PoIChoice,
         estPsi = estPsi,
         scores = scores ,
         estEF = estEF,
         Y_hat = Y_hat, 
         estBeta = beta_hat(grd),
         BIC = LMandBIC$BIC,
         hat_matrix = hat_matrix
         )
}