diff --git a/R/Zalpha_all.R b/R/Zalpha_all.R
index ec55c1c1a6c933c08252e94c91ed05da4004a1b1..35e99918a929496750cfd5fcbb93fec4ef07efb2 100644
--- a/R/Zalpha_all.R
+++ b/R/Zalpha_all.R
@@ -29,7 +29,7 @@
#' @param LDprofile_sd Optional. A numeric vector containing the standard deviation of the \eqn{r^2}{r^2} values for the corresponding bin in the LD profile.
#' @param LDprofile_Beta_a Optional. A numeric vector containing the first estimated Beta parameter for the corresponding bin in the LD profile.
#' @param LDprofile_Beta_b Optional. A numeric vector containing the second estimated Beta parameter for the corresponding bin in the LD profile.
-#' @param minRandL Minimum number of SNPs in each set R and L for the statistics to be calculated. Default is 4.
+#' @param minRandL Minimum number of SNPs in each set R and L for the statistics to be calculated. L is the set of SNPs to the left of the target SNP and R to the right, within the given window size \code{ws}. Default is 4.
#' @param minRL Minimum value for the product of the set sizes for R and L. Default is 25.
#' @param X Optional. Specify a region of the chromosome to calculate the statistics for in the format \code{c(startposition, endposition)}. The start position and the end position should be within the extremes of the positions given in the \code{pos} vector. If not supplied, the function will calculate the statistics for every SNP in the \code{pos} vector.
#'
@@ -213,11 +213,11 @@ Zalpha_all <- function(pos, ws, x=NULL, dist=NULL, LDprofile_bins=NULL, LDprofil
## check L, R and LR
noL <- length(pos[pos>=currentPos-ws/2 & pos < currentPos]) ## Number of SNPs to the left of the current SNP
noR <- length(pos[pos<=currentPos+ws/2 & pos > currentPos]) ## Number of SNPs to the right of the current SNP
+ outputList$LR[i]<-noL*noR
+ outputList$L_plus_R[i]<-choose(noL,2)+choose(noR,2)
if (noL < minRandL || noR < minRandL || noL*noR < minRL){
#NA for everything - leave as is
} else {
- outputList$LR[i]<-noL*noR
- outputList$L_plus_R[i]<-choose(noL,2)+choose(noR,2)
if (is.null(x)==FALSE){
##Left
Lrsq <- lower_triangle(cor(t(x[pos>=currentPos-ws/2 & pos < currentPos,]),use="pairwise.complete.obs")^2)
@@ -276,8 +276,10 @@ Zalpha_all <- function(pos, ws, x=NULL, dist=NULL, LDprofile_bins=NULL, LDprofil
}
}
}
- if (sum(is.na(outputList$LR))==outputLength){
- warning("No statistics were calculated, try reducing minRandL and minRL or increasing the window size")
+ if (length(outputList)>3){
+ if (sum(sapply(outputList[-c(1:3)],function(x) sum(is.na(x))==outputLength))==length(outputList)-3){
+ warning("No statistics were calculated, try reducing minRandL and minRL or increasing the window size")
+ }
}
return(outputList)
}
diff --git a/tests/testthat/test-Zalpha_all.R b/tests/testthat/test-Zalpha_all.R
index 13ec986fb35a76f0811dae6a35ad5ec2bea023c1..8240dbf03a491113c791de35e38940c81005f11c 100644
--- a/tests/testthat/test-Zalpha_all.R
+++ b/tests/testthat/test-Zalpha_all.R
@@ -59,8 +59,8 @@ test_that("Zalpha_all calculates statistics correctly", {
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7]), dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b, minRandL = 4, minRL = 25, X = NULL),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
@@ -83,8 +83,8 @@ test_that("Zalpha_all calculates the statistics correctly with a different windo
expect_equal(Zalpha_all(pos = df$POS, ws = 1100, x = as.matrix(df[,3:7]), dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b, minRandL = 4, minRL = 25, X = NULL),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,NA,25,25,25,25,25,NA,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,NA,20,20,20,20,20,NA,NA,NA,NA,NA),
+ LR=c(0,5,10,15,20,25,25,25,25,25,20,15,10,5,0),
+ L_plus_R=c(10,10,11,13,16,20,20,20,20,20,16,13,11,10,10),
Zalpha_expected=c(NA,NA,NA,NA,NA,0.397641067838740,0.407771198285253,0.412830128247348,0.419512588688985,0.421536156959405,NA,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,NA,0.369069693255773,0.381250676471648,0.385948372025294,0.388272359451172,0.387629142774876,NA,NA,NA,NA,NA),
Zalpha=c(NA,NA,NA,NA,NA,((6+1/4)/10+(2+19/48)/10)/2,((5+5/18)/10+(2+19/48)/10)/2,((2+71/72)/10+(2+19/48)/10)/2,((2+3/16)/10+(2+7/144)/10)/2,((2+3/16)/10+(1+121/144)/10)/2,NA,NA,NA,NA,NA),
@@ -110,8 +110,8 @@ test_that("Zalpha_all calculates all statistics correctly with character matrix"
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df1[,3:7]), dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b, minRandL = 4, minRL = 25, X = NULL),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
@@ -415,8 +415,8 @@ test_that("Zalpha_all calculates statistics correctly with no optional parameter
expect_equal(Zalpha_all(pos = df$POS, ws = 3000),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA)
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91)
))
})
@@ -427,8 +427,8 @@ test_that("Zalpha_all calculates statistics correctly with only x supplied", {
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7])),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
Zbeta=c(NA,NA,NA,NA,(10+5/18)/40,(10+35/36)/45,(11+103/144)/48,(12+73/144)/49,(11+83/144)/48,(11+17/48)/45,(10+65/144)/40,NA,NA,NA,NA)
),tolerance=0.0001)
@@ -441,8 +441,8 @@ test_that("Zalpha_all calculates statistics correctly with dist, LDprofile_bins
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA)
),tolerance=0.0001)
@@ -455,8 +455,8 @@ test_that("Zalpha_all calculates statistics correctly with x, dist, LDprofile_bi
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7]), dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
@@ -475,8 +475,8 @@ test_that("Zalpha_all calculates statistics correctly with only LDprofile_Beta_a
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7]), dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
@@ -497,8 +497,8 @@ test_that("Zalpha_all calculates statistics correctly with only LDprofile_sd not
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7]), dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
@@ -519,8 +519,8 @@ test_that("Zalpha_all calculates statistics correctly with only x not supplied",
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, dist = df$dist, LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha_expected=c(NA,NA,NA,NA,0.390457304338967,0.392014054942343,0.397339546324536,0.398874728980465,0.400715520018796,0.401718327864356,0.399526703832832,NA,NA,NA,NA),
Zbeta_expected=c(NA,NA,NA,NA,0.350404001770323,0.357372168444253,0.360647397851440,0.361951221318345,0.362988750761055,0.364752550557575,0.366343120440209,NA,NA,NA,NA)
),tolerance=0.0001)
@@ -533,8 +533,8 @@ test_that("Zalpha_all calculates statistics correctly with dist not supplied", {
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7]), LDprofile_bins = LDprofile$bin, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
Zbeta=c(NA,NA,NA,NA,(10+5/18)/40,(10+35/36)/45,(11+103/144)/48,(12+73/144)/49,(11+83/144)/48,(11+17/48)/45,(10+65/144)/40,NA,NA,NA,NA)
),tolerance=0.0001)
@@ -547,8 +547,8 @@ test_that("Zalpha_all calculates statistics correctly with LDprofile_bins not su
expect_equal(Zalpha_all(pos = df$POS, ws = 3000, x = as.matrix(df[,3:7]), dist = df$dist, LDprofile_rsq = LDprofile$rsq, LDprofile_sd = LDprofile$sd, LDprofile_Beta_a = LDprofile$Beta_a, LDprofile_Beta_b = LDprofile$Beta_b),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha=c(NA,NA,NA,NA,((3+1/2)/6+(11+41/144)/45)/2,((6+1/4)/10+(9+41/48)/36)/2,((7+31/72)/15+(7+13/48)/28)/2,((8+17/144)/21+(4+7/16)/21)/2,((9+131/144)/28+(2+13/16)/15)/2,((13+97/144)/36+(1+121/144)/10)/2,((15+25/48)/45+(1+55/144)/6)/2,NA,NA,NA,NA),
Zbeta=c(NA,NA,NA,NA,(10+5/18)/40,(10+35/36)/45,(11+103/144)/48,(12+73/144)/49,(11+83/144)/48,(11+17/48)/45,(10+65/144)/40,NA,NA,NA,NA)
),tolerance=0.0001)
@@ -562,8 +562,8 @@ test_that("Zalpha_all calculates statistics correctly with missing value", {
expect_equal(Zalpha_all(pos = df1$POS, ws = 3000, x = as.matrix(df1[,3:7])),
list(
position=c(100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),
- LR=c(NA,NA,NA,NA,40,45,48,49,48,45,40,NA,NA,NA,NA),
- L_plus_R=c(NA,NA,NA,NA,51,46,43,42,43,46,51,NA,NA,NA,NA),
+ LR=c(0,13,24,33,40,45,48,49,48,45,40,33,24,13,0),
+ L_plus_R=c(91,78,67,58,51,46,43,42,43,46,51,58,67,78,91),
Zalpha=c(NA,NA,NA,NA,0.434953703703704,0.473283179012346,0.397114748677249,0.317791005291005,0.300801917989418,0.322897376543210,0.360532407407407,NA,NA,NA,NA),
Zbeta=c(NA,NA,NA,NA,0.248611111111111,0.235185185185185,0.233651620370370,0.257794784580499,0.250144675925926,0.259413580246914,0.271354166666667,NA,NA,NA,NA)
),tolerance=0.0001)