# R/meancoef.R In SpiceFP: Sparse Method to Identify Joint Effects of Functional Predictors

#### Documented in meancoef

#' meancoef
#'
#' This function can be used to compute the mean of coefficients from different
#' partitions in the context of the spicefp approach.
#'
#' @param coef.list list. The second element of the coef_spicefp function
#' outputs. It has the same name as the argument.
#' @param weight a numerical vector of weights with the same length as
#' coef.list.
#'
#' @details  Here, the fine-mesh coefficients are weighted and a weighted mean is deduced.
#' If the user wishes, he can use as weights the slopes associated with the qualities
#' of the models concerned.
#'
#' @return Returns a list of :
#'
#' \describe{
#' \item{weighted_mean}{fine-mesh matrix or array with the weighted mean of the
#'  coefficients}
#' \item{y.estimated}{weighted estimation of \eqn{X \beta}}
#' \item{coefficients.array}{An array with all the fine-mesh coefficients that
#' will be used to compute the weighted mean}
#' \item{weight}{same as inputs}
#' }
#'
#' @importFrom stats weighted.mean
#' @importFrom purrr map_dfc
#'
#' @examples
#' \donttest{
#' ##linbreaks: a function allowing to obtain breaks linearly
#' linbreaks<-function(x,n){
#'      sort(round(seq(trunc(min(x)),
#'                 ceiling(max(x)+0.001),
#'                 length.out =unlist(n)+1),
#'             1)
#'           )
#' }
#' # In this example, we will evaluate 2 candidates with 14 temperature
#' # classes and 15 irradiance classes. The irradiance breaks are obtained
#' # according to a log scale (logbreaks function) with different alpha
#' # parameters for each candidate (0.005, 0.01).
#' ## Data and inputs
#' tpr.nclass=14
#' irdc.nclass=15
#' irdc.alpha=c(0.005, 0.01)
#' p2<-expand.grid(tpr.nclass, irdc.alpha, irdc.nclass)
#' parlist.tpr<-split(p2[,1], seq(nrow(p2)))
#' parlist.irdc<-split(p2[,2:3], seq(nrow(p2)))
#' parlist.irdc<-lapply(
#'    parlist.irdc,function(x){
#'    list(x[],x[])}
#' )
#' m.tpr <- as.matrix(Temperature[,-c(1)])
#'
#' # For the constructed models, only two regularization parameter ratios
#' # penratios=c(1/25,5) is used. In a real case, more candidates
#' # and regularization parameter ratios should be evaluated.

#' ex_sp<-spicefp(y=FerariIndex_Difference$fi_dif, #' fp1=m.irdc, #' fp2=m.tpr, #' fun1=logbreaks, #' fun2=linbreaks, #' parlists=list(parlist.irdc, #' parlist.tpr), #' penratios=c(1/25,5), #' appropriate.df=NULL, #' nknots = 100, #' ncores =2, #' write.external.file = FALSE) #' #' ## Focus on the 2 best models retained by the AIC criterion at iteration 1 #' c.mdls <- coef_spicefp(ex_sp, iter_=1, criterion ="AIC_", #' nmodels=2, ncores = 2, #' dim.finemesh=c(1000,1000), #' write.external.file = FALSE) #' #' # meancoef #' # Compute the mean of the coefficients of these models #' mean.c.mdls<-meancoef(c.mdls$coef.list,
#'                      weight = c.mdls$Model.parameters$Slope_)
#' g3<-mean.c.mdls$weighted_mean #' g3.x<-as.numeric(rownames(g3)) #' g3.y<-as.numeric(colnames(g3)) #' #' #' library(fields) #' plot(c(10,2000),c(15,45),type= "n", axes = FALSE, #' xlab = "Irradiance (mmol/m2/s - Logarithmic scale)", #' ylab = "Temperature (deg C)",log = "x") #' rect(min(g3.x),min(g3.y),max(g3.x),max(g3.y), col="black", border=NA) #' image.plot(g3.x,g3.y,g3, horizontal = FALSE, #' col=designer.colors(256, c("blue","white","red")), #' add = TRUE) #' axis(1) ; axis(2) #' #' closeAllConnections() #' #'} #' #' @export meancoef<-function(coef.list, weight){ try(if(length(weight)!=length(coef.list)) stop("The length of weight vector must be equal to the length of coef.list")) try(if(sum(is.na(weight))!=0) stop("weight vector must not contain missing data")) dim.finemesh=dim(coef.list[]$Candidate.coef.NA.finemeshed)
coef.NA.finemeshed=list()
for (g in 1:length(coef.list)){
coef.NA.finemeshed[[g]]<-coef.list[[g]]$Candidate.coef.NA.finemeshed } spicefp_coefweimean<-function(X,W){ if(setequal(X, rep(NA,length(X)))) {NA} else { X[which(is.na(X))]<-0 weighted.mean(X,W) } } if(length(dim(coef.list[]$Candidate.coef.NA.finemeshed))==2){
acoef<-array(NA,c(length(coef.NA.finemeshed),
dim.finemesh,
dim.finemesh))
for (i in 1:length(coef.NA.finemeshed)) {
acoef[i,,]<-coef.NA.finemeshed[[i]]
}
spicefp.coef= apply(acoef,
c(2,3),
FUN = function(x){spicefp_coefweimean(X=x, W=weight)})
dimnames(spicefp.coef)=dimnames(coef.NA.finemeshed[])
} else {
acoef<-array(NA,c(length(coef.NA.finemeshed),
dim.finemesh,
dim.finemesh,
dim.finemesh))
for (i in 1:length(coef.NA.finemeshed)) {
acoef[i,,,]<-coef.NA.finemeshed[[i]]
}
spicefp.coef= apply(acoef,
c(2,3,4),
FUN = function(x){spicefp_coefweimean(X=x, W=weight)})
dimnames(spicefp.coef)=dimnames(coef.NA.finemeshed[])
}

suppressMessages(
y_est <-apply(map_dfc(1:length(coef.list),
function(x){coef.list[[x]]\$y.estimated*weight[x]}),1,
sum)/sum(weight)
)

return(list(weighted_mean=spicefp.coef,
y.estimated=y_est,
coefficients.array=acoef,
weight=weight))
}


## Try the SpiceFP package in your browser

Any scripts or data that you put into this service are public.

SpiceFP documentation built on Sept. 15, 2021, 9:07 a.m.