Merge pull request #1 from DecodeGenetics/feature/fixing_checks

No issues in win dev or rhub

Author: gumeo

.Rbuildignore

@@ -1,3 +1,5 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 ^LICENSE\.md$
+^README.Rmd
+^cran-comments\.md$

DESCRIPTION

@@ -1,21 +1,24 @@
 Package: gnonadd
 Type: Package
 Title: Various Non-Additive Models for Genetic Associations
-Version: 1.0.0
+Version: 1.0.1
 Authors@R: c(
      person("Audunn S.", "Snaebjarnarson", , "audunn.snaebjarnarson@decode.is", role = c("aut", "cre")),
      person("Gudmundur", "Einarsson", , "gudmundur.einarsson2@decode.is", role = c("aut")),
      person("Daniel F.", "Gudbjartsson", , "daniel.gudbjartsson@decode.is", role = c("aut"))
      )
 Description: The goal of gnonadd is to simplify workflows in the analysis of non-additive effects of 
     sequence variants. This includes variance effects, correlation effects, interaction 
-    effects and dominance effects.
+    effects and dominance effects. The package also includes convenience functions for visualization.
+URL: https://github.com/DecodeGenetics/gnonadd
+BugReports: https://github.com/DecodeGenetics/gnonadd/issues
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.2.1
 Depends: 
     R (>= 2.10)
 Imports: 
-    ggplot2,
+    ggplot2
+Suggests:
     MASS

NEWS.md

@@ -0,0 +1,4 @@
+# gnonadd 1.0.1
+
+* Added a `NEWS.md` file to track changes to the package.
+* Making package ready for initial CRAN release

R/Dominance_CaseControl_model.R

@@ -13,7 +13,7 @@
 #' @param sex A numerical vector containing sex, coded 0 for males, 1 for females and -1 for unknown
 #' @param round_imputed A boolian variable determining whether imputed genotype values should be rounded to the nearest integer in the analysis
 #' @param covariates A dataframe containing any other covariates that should be used; one column per covariate.
-#' 
+#'
 #' @returns
 #' A list with the dominanc effect (on log-scale) and corresponding standard error, z statistic and p-value
 #' @examples
@@ -39,11 +39,11 @@ dominance_CC.calc <- function(cc, g, yob=rep(-1,length(cc)), sex=rep(-1,length(c
     sex <- as.factor(sex)
     g <- g - mean(g)
     g2 <- as.numeric(g_rounded == 2)
-    
+
     #We define a dataframe containing all variables that should be considered
     Dom_data <- as.data.frame(cbind(cc, g2))
     Dom_data <- cbind(Dom_data, g)
-    if(sd(yob) > 0) {
+    if(stats::sd(yob) > 0) {
       Dom_data <- cbind(Dom_data, yob)
     }
     if(length(unique(no_date)) > 1) {
@@ -55,15 +55,15 @@ dominance_CC.calc <- function(cc, g, yob=rep(-1,length(cc)), sex=rep(-1,length(c
     if(nrow(covariates) > 0) {
       Dom_data <- cbind(Dom_data, covariates)
     }
-    
+
     #We use logistic regression to estimate the dominance effect
-    l_delta <- glm(cc ~ ., data = Dom_data, family = 'binomial')
+    l_delta <- stats::glm(cc ~ ., data = Dom_data, family = 'binomial')
     param <- "g2"
-    if(param %in% rownames(coef(summary(l_delta)))){
+    if(param %in% rownames(stats::coef(summary(l_delta)))){
       delta <- summary(l_delta)$coeff[param, 1]
       se <- summary(l_delta)$coeff[param, 2]
       z <- summary(l_delta)$coeff[param, 3]
-      p <- summary(l_delta)$coeff[param, 4] 
+      p <- summary(l_delta)$coeff[param, 4]
     }else{
       warning("Singular model matrix")
       delta <- NA

R/Dominance_model.R

@@ -11,7 +11,7 @@
 #' @param g A vector with (possibly imputed) genotype values. All entries should be larger than 0 and smaller than 2.
 #' @param round_imputed A boolian variable determining whether imputed genotype values should be rounded to the nearest integer in the analysis
 #' @param covariates A dataframe containing any covariates that should be used; one column per covariate.
-#' 
+#'
 #' @returns
 #' A list with the dominanc effect and corresponding standard error, t statistic and p-value
 #' @examples
@@ -32,22 +32,22 @@ dominance.calc <- function(qt, g, round_imputed = F, covariates = as.data.frame(
     p <- NA
   }else {
     g2 <- as.numeric(g_rounded == 2)
-    
+
     #We define a dataframe containing all variables that should be considered
     Dom_data <- as.data.frame(cbind(qt, g2))
     Dom_data <- cbind(Dom_data, g)
     if(nrow(covariates) > 0) {
       Dom_data <- cbind(Dom_data, covariates)
     }
-    
+
     #We use linear regression to estimate the dominance effect
-    l_delta <- lm(qt ~ ., data = Dom_data)
+    l_delta <- stats::lm(qt ~ ., data = Dom_data)
     param <- "g2"
-    if(param %in% rownames(coef(summary(l_delta)))){
+    if(param %in% rownames(stats::coef(summary(l_delta)))){
       delta <- summary(l_delta)$coeff[param, 1]
       se <- summary(l_delta)$coeff[param, 2]
       t <- summary(l_delta)$coeff[param, 3]
-      p <- summary(l_delta)$coeff[param, 4] 
+      p <- summary(l_delta)$coeff[param, 4]
     }else{
       warning("Singular model matrix")
       delta <- NA

R/Ellipse_by_genotype.R

@@ -24,6 +24,13 @@
 #' @export ellipse.by.gen
 ellipse.by.gen <- function(qt1, qt2, g, trait_name1 = 'qt trait 1', trait_name2 = 'qt trait 2',
                         title = '', sample_size = 500) {
+  g_factor <- NULL
+  x0 <- NULL
+  x1 <- NULL
+  x2 <- NULL
+  y0 <- NULL
+  y1 <- NULL
+  y2 <- NULL
   g <- round(g)
   D <- cbind(qt1, qt2)
   D <- cbind(D, g)
@@ -42,7 +49,7 @@ ellipse.by.gen <- function(qt1, qt2, g, trait_name1 = 'qt trait 1', trait_name2
       D_sample <- rbind(D_sample, D_temp[sample(1:nrow(D_temp), size = min(sample_size, nrow(D_temp)), replace = FALSE), ])
       qt1_mean <- mean(D_temp$qt1)
       qt2_mean <- mean(D_temp$qt2)
-      Sigma <- cov(D_temp[, c(1, 2)])
+      Sigma <- stats::cov(D_temp[, c(1, 2)])
       Princip <- eigen(Sigma)
       flip_direction1 <- 0
       flip_direction2 <- 0
@@ -66,7 +73,7 @@ ellipse.by.gen <- function(qt1, qt2, g, trait_name1 = 'qt trait 1', trait_name2
   }
   ggplot2::ggplot(D_sample, ggplot2::aes(x = qt1 , y = qt2 ,color = g_factor))+
     ggplot2::geom_point()+ggplot2::theme_classic()+
-    ggplot2::geom_smooth(method = 'lm', data = D, se = F, formula = as.formula('y ~ x')) +
+    ggplot2::geom_smooth(method = 'lm', data = D, se = F, formula = stats::as.formula('y ~ x')) +
     ggplot2::coord_fixed() +
     ggplot2::scale_color_manual(values = c('Non-carriers' = '#F8766D', 'Heterozygotes' = '#00BA38', 'Homozygotes' = '#619CFF')) +
     ggplot2::geom_segment(ggplot2::aes(x = Arrow_data[1, 1], y = Arrow_data[1, 2],

R/Env_interaction_CaseControl_all_vs_all.R

@@ -2,7 +2,7 @@
 #'
 #' @description
 #' Given a set of variants and environmental traits, and a single case control variable, this function calculates the interaction effect of all possible variant-environmental pairs
-#' 
+#'
 #' @param cc A numeric vector
 #' @param g A matrix, where each colomn represents a variant
 #' @param env A matrix, where each row represents an environmental variable
@@ -18,16 +18,19 @@
 #' @returns
 #' A dataframe with all possible variant-environmental pairs and their estimated interaction effect
 #' @examples
-#' g_vec <- matrix(0, nrow = 100000, ncol = 3)
+#' N_run <- 25000
+#' g_vec <- matrix(0, nrow = N_run, ncol = 3)
 #' freqs <- runif(ncol(g_vec), min = 0, max = 1)
-#' env_vec <- matrix(0, nrow = 100000, ncol = 3)
+#' env_vec <- matrix(0, nrow = N_run, ncol = 3)
 #' for(i in 1:ncol(g_vec)){
-#'  g_vec[, i] <- rbinom(100000, 2, freqs[i])
+#'  g_vec[, i] <- rbinom(N_run, 2, freqs[i])
 #' }
 #' for( i in 1:ncol(env_vec)){
-#'  env_vec[, i] <- round(runif(100000,min=0,max=6))
+#'  env_vec[, i] <- round(runif(N_run,min=0,max=6))
 #' }
-#' cc_vec <- rbinom(100000,1,0.1 * (1.05 ^ g_vec[, 1]) * (1.06 ^ env_vec[,1]) * (0.95 ^ g_vec[, 2]) * (1.1^(g_vec[, 1] * env_vec[, 1])))
+#' cc_vec <- rbinom(N_run,1,0.1 * (1.05 ^ g_vec[, 1]) *
+#'           (1.06 ^ env_vec[,1]) * (0.95 ^ g_vec[, 2]) *
+#'           (1.1^(g_vec[, 1] * env_vec[, 1])))
 #' res <- pairwise_env_int_CC.calc(cc_vec, g_vec, env_vec)
 #' @export
 pairwise_env_int_CC.calc <- function(cc, g, env, yob = rep(-1,length(cc)), sex = rep(-1,length(cc)),
@@ -44,7 +47,7 @@ pairwise_env_int_CC.calc <- function(cc, g, env, yob = rep(-1,length(cc)), sex =
       counter <- counter + 1
       A$variant_name[counter] <- variant_names[i]
       A$env_name[counter] <- env_names[j]
-      res <- env_interaction_CC.calc(cc, g[, i], env[, j], yob = yob, sex = sex, 
+      res <- env_interaction_CC.calc(cc, g[, i], env[, j], yob = yob, sex = sex,
                                      round_imputed = round_imputed, dominance_term = dominance_term, square_env = square_env, covariates = covariates )
       A$int_effect[counter] <- res$interaction_effect
       A$se[counter] <- res$standard_error

R/Env_interaction_model.R

@@ -16,7 +16,7 @@
 #' @param dominance_term A boolian variable determining whether a dominance term for the variant should be included as a covariates in the analysis
 #' @param square_env A boolian variable determining whether the square of the environmental trait should be included as a covariate in the analysis
 #' @param covariates A dataframe containing any other covariates that should be used; one column per covariate
-#' 
+#'
 #' @returns
 #' A list with the environmental interaction effect and corresponding standard error, t statistic and p-value
 #' @examples
@@ -28,19 +28,19 @@
 env_interaction.calc <- function(qt, g, env, round_imputed = F, dominance_term = F,
                                  square_env = F,  covariates = as.data.frame(matrix(0, nrow = 0, ncol = 0))){
   r <- rank(env)
-  env_normal <- qnorm(r / (length(r) + 1))
+  env_normal <- stats::qnorm(r / (length(r) + 1))
   if(round_imputed == T){
     g <- round(g)
   }
   int <- g * env_normal
-  if(sd(int) == 0){
+  if(stats::sd(int) == 0){
     warning("Interaction undefined. All interaction values are the same.")
     gamma <- NA
     se <- NA
     t <- NA
     p <- NA
   }else{
-    
+
     #We define a dataframe containing all variables that should be considered
     Env_int_data <- as.data.frame(cbind(qt, int))
     Env_int_data <- cbind(Env_int_data, g)
@@ -54,15 +54,15 @@ env_interaction.calc <- function(qt, g, env, round_imputed = F, dominance_term =
     if(nrow(covariates) > 0) {
       Env_int_data <- cbind(Env_int_data, covariates)
     }
-    
+
     #We use linear regression to estimate the environmental interaction effect
-    l_interaction <- lm(qt ~ ., data = Env_int_data)
+    l_interaction <- stats::lm(qt ~ ., data = Env_int_data)
     param <- "int"
-    if(param %in% rownames(coef(summary(l_interaction)))){
+    if(param %in% rownames(stats::coef(summary(l_interaction)))){
       gamma <- summary(l_interaction)$coeff[param, 1]
       se <- summary(l_interaction)$coeff[param, 2]
       t <- summary(l_interaction)$coeff[param, 3]
-      p <- summary(l_interaction)$coeff[param, 4] 
+      p <- summary(l_interaction)$coeff[param, 4]
     }else{
       warning("Singular model matrix")
       gamma <- NA

R/Env_interactions_all_vs_all.R

@@ -27,7 +27,9 @@
 #'  env_vec[, i] <- round(runif(100000,min=0,max=6))
 #' }
 #'
-#' qt_vec <- rnorm(100000) + 0.1 * g_vec[, 1] + 0.2 * g_vec[, 2] -0.1 * env_vec[, 3] + 0.1 * env_vec[, 1] + 0.1 * g_vec[, 1] * env_vec[, 1]
+#' qt_vec <- rnorm(100000) + 0.1 * g_vec[, 1] + 0.2 *
+#'           g_vec[, 2] -0.1 * env_vec[, 3] + 0.1 *
+#'           env_vec[, 1] + 0.1 * g_vec[, 1] * env_vec[, 1]
 #' res <- pairwise_env_int.calc(qt_vec, g_vec, env_vec)
 #' @export
 pairwise_env_int.calc <- function(qt, g, env, round_imputed = F, dominance_term = F, square_env = F,

R/Histogram_by_genotype.R

@@ -21,16 +21,17 @@
 #' hist_by_gen(qt_vec, geno_vec)
 #' @export hist_by_gen
 hist_by_gen <- function(qt, g, bins = 100, trait_name = 'qt trait', title = '', outlier_quantiles = c(0.025, 0.975), sd_lines = c(1,1)) {
+  outlier <- NULL
   g <- round(g)
   D <- as.data.frame(cbind(qt,g))
   D$g_factor <- factor(D$g,levels=0:2, labels = c('Non-carriers', 'Heterozygotes', 'Homozygotes'))
-  quant_vals <- as.numeric(quantile(D$qt, outlier_quantiles))
+  quant_vals <- as.numeric(stats::quantile(D$qt, outlier_quantiles))
   D$outlier <- as.factor(0+(D$qt<quant_vals[1] | D$qt > quant_vals[2]))
   a <- sd_lines[1]
   b <- sd_lines[2]
   vertical_lines <- data.frame(g_factor = levels(D$g_factor), mean = c(mean(D$qt[D$g == 0]), mean(D$qt[D$g == 1]), mean(D$qt[D$g == 2])),
-                         lower_sd <- c(mean(D$qt[D$g == 0]) - a * sd(D$qt[D$g == 0]), mean(D$qt[D$g == 1]) - a * sd(D$qt[D$g == 1]), mean(D$qt[D$g == 2]) - a * sd(D$qt[D$g == 2])),
-                         upper_sd <- c(mean(D$qt[D$g == 0]) + b * sd(D$qt[D$g == 0]), mean(D$qt[D$g == 1]) + b * sd(D$qt[D$g == 1]), mean(D$qt[D$g == 2]) + b * sd(D$qt[D$g == 2])))
+                         lower_sd <- c(mean(D$qt[D$g == 0]) - a * stats::sd(D$qt[D$g == 0]), mean(D$qt[D$g == 1]) - a * stats::sd(D$qt[D$g == 1]), mean(D$qt[D$g == 2]) - a * stats::sd(D$qt[D$g == 2])),
+                         upper_sd <- c(mean(D$qt[D$g == 0]) + b * stats::sd(D$qt[D$g == 0]), mean(D$qt[D$g == 1]) + b * stats::sd(D$qt[D$g == 1]), mean(D$qt[D$g == 2]) + b * stats::sd(D$qt[D$g == 2])))
 
   ggplot2::ggplot(D, ggplot2::aes(x = qt, fill = outlier)) + ggplot2::geom_histogram(color = 'black', bins = bins) + ggplot2::theme_classic() +
     ggplot2::facet_grid(g_factor ~ . , scales = "free_y") + ggplot2::xlab(trait_name) + ggplot2::ggtitle(title) +