### wrapper for MCMCglmm::rbv that allows 0 variances
rbv0 <- (pedigree, G){
X <- (0, =(pedigree), =(G))
index <- ((G)!=0)
((index)>0){
(((G)==0)) G <- G[index,index]
X2 <- MCMCglmm::rbv(pedigree=pedigree, G=G)
X,index] <- X2
}
X
}
### wrapper for mvnfast::rmvn that allows 0 variances
rmvn0 <- (n,mu,){
X <- (0, =n, =())
index <- (()!=0)
X,index] <- mvnfast::rmvn(n=n, mu=mu[index], =[index,index])
X
}
### function to generate ar1 matrix
ar1_cor <- (n, rho) {
exponent <- ((1:n - 1, = n, = n, byrow = ) -
(1:n - 1))
rho^exponent
}
# from https://www.r-bloggers.com/2020/02/generating-correlation-matrix-for-ar1-model/
#ar1_cor(n=10,rho=0.5)
## cov_str_check
## performs error checking of all cov structures
cov_str_check <- (data_structure, pedigree, phylogeny, cov_str, parameters,){
cs_check <- (("pedigree","phylogeny","cov_str"), (j){
cs <- (j)
(!(cs) | (cs)) (j, " needs to be a list", call.=)
})
param_names <- (parameters)
group_names <- (parameters,(x) x$)
ped_names <- (pedigree)
phylo_names <- (phylogeny)
cov_names <- (cov_str)
cs_names <- (ped_names,phylo_names,cov_names)
(((cs_names))) ("Cannot have multiple covariance structures (pedigree/phylogeny/cov_str) linking to the same item in the parameter list. If multiple covariance structures are needed to be linked to the same grouping factor in the data_structure (for example simulating additive genetic and dominance effects), then create multiple items in the parameter list, with different names, but the same 'group', and link the covariance structures accordingly.", call.=)
((!cs_names %in% param_names)) ("Some names in pedigree/phylogeny/cov_str are not in the parameter list", call.=)
((data_structure), (i){
# i = colnames(data_structure)[1]
# are any of the parameter list names associated with it
list_names <- param_names[group_names %in% i]
ped_link <- list_names[list_names %in% ped_names]
phylo_link <- list_names[list_names %in% phylo_names]
cov_link <- list_names[list_names %in% cov_names]
all_link <- (ped_link,phylo_link,cov_link)
((all_link)>1){
("Multiple covariance structures linked to ",i,". The function assumes that the covariance structures are ordered exactly the same. If they are not then the simulations will not run as you expect. You will need to create multiple columns in grouping structure and link different covariance structures to each one.", call.=)
}
})
# ped_names <- c(names(pedigree),names(phylogeny),names(cov_str))
# names_check <- lapply(ped_names,function(i){
# # data_structure[,i]
# if(!all(unique(rownames(chol_str_all[[i]])) %in% unique(data_structure[,parameters[[i]]$group]))) stop(paste("all IDs in the pedigree/phylogeny/cov_str linked with", i, "are not in the data_structure"), call.=FALSE)
# if(!all(unique(data_structure[,parameters[[i]]$group]) %in% unique(rownames(chol_str_all[[i]])))) stop(paste("all IDs in data_structure are not in the pedigree/phylogeny/cov_str linked with", i), call.=FALSE)
# })
}
### function to turn data_structure into indexes. Matches names in data_structure to linked pedigree/phylogeny/cov_str to make sure indexing is correct. Doesnt do any error checking
index_factors <- (data_structure, pedigree, phylogeny, cov_str, parameters, index_link,suppress_index_warning,){
p_names <- (parameters)!(parameters)%in%("intercept","interactions")]
group_names <- (p_names,(x) parameters[x]]$)
new_ds <- ((data_structure)){
} ((pedigree) & (phylogeny) & (cov_str)){
(data_structure,2,(x) ((x)))
}{
(,
## for each column of the data structure
((data_structure), (i){
# i = colnames(data_structure)[1]
# are any of the parameter list names associated with it
list_names <- group_names[group_names %in% i]
ped_link <- list_names[list_names %in% (pedigree)]
phylo_link <- list_names[list_names %in% (phylogeny)]
cov_link <- list_names[list_names %in% (cov_str)]
# if linked with a something get the row number in the relevant cov_str, so the indexing will match the order in that cov str.
# Can only match with one so is assuming that if multiple things are linked with it that the ordering is the same, so takes the first linked cov str it can find and indexes according to that
((ped_link)>0){
# first column of pedigree
(data_structure,i],pedigree[ped_link[1]]],1])
} ((phylo_link)>0){
# names of phylogeny
(data_structure,i],phylogeny[phylo_link[1]]]$tip.label)
} ((cov_link)>0){
# rownames(of cov_str)
(data_structure,i],(cov_str[cov_link[1]]]))
}{
((data_structure,i]))
}
})
)
}
(new_ds) <- (data_structure)
(!(index_link)){
more_ds <- (,(index_link, (x){
ds_links <- (x, '-')[1]]
## give warning if not all levels match
(!(data_structure,ds_links[1]] %in% data_structure,ds_links[2]]) & !suppress_index_warning) (("Not all levels are of", ds_links[1], "are present in", ds_links[2], "meaning that there will be NAs in the new grouping factor"), call.=)
new_ds,ds_links[2]][match(data_structure,ds_links[1]], data_structure,ds_links[2]])]
}
))
((new_ds,more_ds))
}{
(new_ds)
}
}
cov_str_list <- (parameters, data_structure, phylogeny, phylogeny_type, cov_str,){
#pedigree, pedigree_type,
# ped_chol <- sapply(names(pedigree), function(x){
# if(pedigree_type[[x]]=="A") Matrix::chol(nadiv::makeA(pedigree[[x]]))
# else if(pedigree_type[[x]]=="D") Matrix::chol(nadiv::makeD(pedigree[[x]]))
# else if(pedigree_type[[x]]=="E") Matrix::chol(nadiv::makeAA(pedigree[[x]]))
# })
phylo_chol <- ((phylogeny), (x){
phylo_vcv <- ape::vcv(phylogeny[x]], corr = )
# if(phylogeny_type[[x]]=="OU") {
# ## need way of specifying alpha
# phylo_vcv <- exp(phylo_vcv * - alpha)
# diag(phylo_vcv) <- 1
# }
::((phylo_vcv), "dgCMatrix")
})
cor_chol <- (cov_str, (x) ::((x), "dgCMatrix"))
chol_str_all<-(phylo_chol,cor_chol)#ped_chol,
add_list<-(parameters)!(parameters) %in% ((chol_str_all),"intercept","interactions")]
(i add_list){
chol_str_all[i]] <- #Matrix::Diagonal(parameters[[i]][["n_level"]])
}
( chol_str_all)
}
sim_predictors <- (parameters,data_structure, str_index, cov_str_all, known_predictors, pedigree, ){
traits <- (, ( (parameters)!(parameters)%in%("intercept","interactions")], (i){
# i<-"animal"
p <- parameters[i]]
k <- (p$)
n <- p$n_level
(p$fixed){
## if factor make design matrix
# x<-stats::model.matrix(stats::formula(paste("~ factor(",p$group,")-1")),as.data.frame(str_index))
## if there are group names in the data structure, then use factor levels with the same order as the specified names
ds_levels_group <- (data_structure,p"group"]]])
((ds_levels_group %in% 1:(ds_levels_group)) ){
fac_levels <- 1:(ds_levels_group)
}{
fac_levels <- p$
}
x<-::(::(("~ factor(",p$,",levels=fac_levels)-1")),(data_structure))
} (p$covariate){
## if covariate, make design matrix from data structure
x<- ((str_index,p$],k),(str_index),k)
}{
## otherwise simulate 'traits' at each level from multivariate normal
# x <- if(is.null(cov_str_all[[i]])) {
# mvnfast::rmvn(n=n, mu=p$mean, sigma=p$vcov)
# }else{
# methods::as(Matrix::crossprod(cov_str_all[[i]],mvnfast::rmvn(n=n, mu=p$mean, sigma=p$vcov)),"matrix")
# }
(i %in% (pedigree)){
## if name is listed in pedigree argument, link to pedigree
x <- rbv0(pedigree[i]],p$)
}{
x <- rmvn0(n=n, mu=p$, =p$)
}
(!(cov_str_all[i]])) {
x <- ::(Matrix::(cov_str_all[i]],x),"matrix")
}
### apply functions
### add them in likes betas, making sure they are in the right order? then apply them to the cols?
x <- (1:k,(i) (x,i],p"functions"]][i]))
## expand traits to be the same length as the number of observations using data structure
(!p$ %in% ("observation","residual")) x <- x[str_index,p$],,=]
}
## use names form parameter list
(x) <- p$
(x)
}))
# traits <- data.frame(a=rnorm(100),b=rnorm(100),c=rnorm(100),d=rnorm(100))
# interactions <- list(names=c("a:b","c:d", "exp(a)"), beta=c(2,3))
## add in known predictors
(!(known_predictors)){
traits<-(traits,known_predictors$predictors)
}
## add in interactions
(!(parameters"interactions"]])){
x_int <- (,((parameters"interactions"]]$,":"), (j){
((=(j, collapse="*")), envir = (traits) )
}))
(x_int) <- parameters"interactions"]]$
traits<-(traits,x_int)
}
(traits)
}
generate_y <- (predictors, intercepts, betas, str_index, model, y_pred_names,extra_param,){
## evaluate model
## - if model is missing, add all simulated predictors together
((model)) {
y <- predictors %*% betas
y <- ((y) + intercepts)
} {
## for evaluation with model formula
beta_predictors <- predictors * (betas, ((predictors),(betas))) #predictors %*% diag(as.vector(betas))
y_predictors <- (beta_predictors,predictors,str_index)
(y_predictors) <- y_pred_names
## allow I() and subsets to be properly linked to y_predictors
model <- ("I\\((\\w+)\\)","\\1_raw",model)
model <- ("\\[(\\w+)\\]","\\[\\1_ID\\]",model)
# evaluate the formula in the context of y_predictors and the extra params
# y <- eval(parse(text=model), envir = c(as.data.frame(y_predictors),as.list(extra_param)))
model2 <- (model,";\n return(data.frame(mget(ls()[!ls() %in% c(colnames(y_predictors),names(extra_param),'intercept')])))")
y <- ((=model2), envir = ((y_predictors),intercept=intercepts,(extra_param)))
(!("intercept",model)) y <- y + intercepts
((y)) y <- (y)
}
(y)
}
generate_y_list <- (parameters, str_index, predictors, model,known_predictors,){
((str_index)) str_index <-
intercepts <- parameters"intercept"]]
parameters <- parameters!(parameters) %in% ("intercept")]
## put all betas together
#order betas to match predictors
betas <- (,(parameters,(x) x$))
(!(known_predictors)){
betas<-(betas,known_predictors$)
}
betas<-betas[colnames(predictors[1]]),]
y_pred_names <- ((predictors[1]]), ((predictors[1]]),"_raw"), (!(str_index)){((str_index),"_ID")})
## extract and name extra parameters
(!(model)){
param_names <- ("names", "group", "mean", "vcov", "vcorr", "beta", "fixed", "covariate", "n_level", "functions")
extra_param <- ((parameters, (x){ x!(x) %in% param_names] }))
(!(extra_param)){
(extra_param) <- ((parameters, (x) (x)!(x) %in% param_names]
))
## check extra param names dont clash with y_trait names
(((extra_param) %in% (y_pred_names))) ("You cannot name extra parameters the same as any variables")
}
}
y <- (predictors, (x) generate_y(x, intercepts=intercepts, betas=betas, str_index=str_index, model=model, y_pred_names=y_pred_names,extra_param=extra_param))
(y)
}
transform_dist <- (y, , link, response_names,){
inv <- (x) 1/x
j <- (y)
((link)==1 & j>1) link <- (link,j)
(()==1 & j>1) <- (,j)
## convert the link argument into an actual function
link_function <-
(link=="log", "exp",
(link=="inverse", "inv",
ifelse (link=="logit", "plogis",
ifelse (link=="probit", "pnorm",
link))))
y_family <- (1:j,(i){
## apply link function to y
y_link <- (link_function[i])(y,i])
## sample from poisson or binomial
([i]=="gaussian") y_link
([i]=="poisson") ::((y_link),y_link)
([i]=="binomial") ::((y_link),1,y_link)
})
## think colnames(y) will be null unless model is specified but should check :)
(((y))){
(!(response_names)){
(y_family) <- response_names
}{
(y_family) <- (j==1)"y" ("y",1:j)
}
}{
(y_family) <- (y)
}
(y_family)
}
