You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Parameter vector to be estimated (first entry is intercept)
beta<- c(0, 1, 1, -1, 0, 0, 1, -1, -1, 0, 0) /4
Creating a simulation model with randomly selected families
make_sim_model<-function(covariates, family_set) {
n<- nrow(covariates)
# 5-dimensional R-vine structure matrixd<-5Matrix<- c(5, 2, 4, 1, 3,
0, 2, 4, 3, 1,
0, 0, 1, 4, 3,
0, 0, 0, 4, 3,
0, 0, 0, 0, 3)
Matrix<-matrix(Matrix, d, d)
# family matrixfamily<-matrix(0, d, d)
# initialize modelcount<-1model<- vector(mode="list", length=d* (d-1) /2)
# model formula for true modeltmpform<-"B.1 + B.2 + B.3 + Z.1 + Z.2 + Z.3"for (iin1:3) {
tmpform<- paste0(tmpform, " + s(S.", i, ", k = 10, bs = 'cr')")
}
# select a family for each pair copula and set GAM relationship to true modelfor (treein1:(d-1)) {
for (pcin1:(d-tree)) {
# select a copula familyfamily[d-tree+1, pc] <- sample(family_set, 1)
# spline approximation of the true functionsS<-covariates[, substr(colnames(covariates), 1, 1) =="S"]
y<- rowSums(sapply(1:5, function(i) s_funs[[i]](S[, i])))
form<- as.formula(paste0("y ~", tmpform))
approx<-mgcv::gam(form, data= cbind(y, covariates))
# create a dummy gamBiCop objecttmp<- gamBiCopFit(
data= cbind(u1= runif(n), u2= runif(n), covariates),
formula=form,
family=1,
n.iters=1
)$res# set to true model
attr(tmp, "model") <-approx
attr(tmp, "family") <-family[d-tree+1, pc]
tmp@model$coefficients[1] <-0tmp@model$coefficients[2:7] <-beta[beta!=0]
if (family[d-tree+1, pc] ==2) {
attr(tmp, "par2") <-4
}
model[[count]] <-tmpcount<-count+1
}
}
# return gamVine model
gamVine(
Matrix=Matrix,
model=model,
covariates= colnames(covariates)
)
}
Extracting results from fitted models
# evaluation gridgrid<- cbind(
matrix(1, 50, 10),
sapply(1:5, function(i) seq(0, 1, l=50))
)
colnames(grid) <- c(paste0("B.", 1:5), paste0("Z.", 1:5), paste0("S.", 1:5))
# extract results for a single pair-copula fitget_pc_results<-function(x) {
gamBiCopPredict(x, newdata=grid, type="terms", target="calib") %>%
map_df(~ tidy_terms(.)) %>%
mutate(family=if (is.list(x)) x$familyelsex@family)
}
# tidy results returned from gamBiCopPredicttidy_terms<-function(x) {
# rename smooths to bare variable names
colnames(x) <- gsub("(^s\\()|(\\))", "", colnames(x))
# add intercept as columnx<- as_tibble(cbind(intercept= attr(x, "constant"), x))
# fill missing columns with 0x[setdiff(colnames(grid), colnames(x))] <-0# standardize variable orderx<-x[c("intercept", colnames(grid))]
x %>%
# in case one of the smooth covariates is selected as a linear term, # it will not integrate to zero; hence, we need to subsume the # means of all S.i into the intercept, and standardize the variables# accordingly
mutate(intercept=intercept- mean(S.1+S.2+S.3+S.4+S.5)) %>%
mutate_at(paste0("S.", 1:5), funs(.- mean(.))) %>%
# add grid sequence for smooth covariates
mutate(t=grid[, 15])
}
res_linear<-res %>%
# estimates for linear coefficients are constant for all t, we can focus on # first value
filter(t==0) %>%
# select linear covariates
select(method, n, tree, seed, starts_with("B"), starts_with("Z")) %>%
gather(term, value, -method, -n, -tree, -seed) %>%
# rename terms to appropriate index of beta
mutate(j= (substr(term, 1, 1) =="Z") *5+ as.numeric(substr(term, 3, 3))) %>%
# calculate mean estimates and quantiles
group_by(j, n, method, tree) %>%
summarize(
mu= mean(value),
q05= quantile(value, 0.05),
q95= quantile(value, 0.95),
true=beta[1+ first(j)]
) %>%
ungroup()