Error in `smooth.spline(lambda, pi0, df = smooth.df)`: missing or infinite values in inputs are not allowed

m <- 300 # this has to be at least this large or `pip` (really `qvalue::pi0est`) fails in stupid ways
# always remove zero variance rows, draw them again
# NOTE: failure to do this results in all kinds of testing errors, like:
# - FSTAT: essentially perfect fit: summary may be unreliable
# - jackstraw_subspace/svd: Error: infinite or missing values in 'x'
# - jackstraw_rpca/rsvd: Error: NA/NaN/Inf in foreign function call (arg 1)
# - jackstraw_irlba/irlba: Error: missing value where TRUE/FALSE needed
# - jackstraw_kmeans/kmeans: Error: NA/NaN/Inf in foreign function call (arg 1)
# - jackstraw_kmeanspp/ClusterR::KMeans_rcpp: Error: the data includes NaN's or +/- Inf values
# - jackstraw_MiniBatchKmeans/ClusterR::MiniBatchKmeans: Error: the data includes NaN's or +/- Inf values
# - jackstraw_pam/pam: Error: No clustering performed, NAs in the computed dissimilarity matrix.
# for PCA and other analyses for continuous data, makes sense to centerscale data
## # since LFA is imported, might as well use it to get actual LFs for our random data
# to have this work without Bioconductor, make fake LFs using PCA, in the same format as LFA's
LF1 <- cbind( eigen( crossprod( Xc ) )$vectors[ , 1:(d-1) ], 1 )
# in practice LF0 is just intercept model
# set these parameters (match defaults), which also determine dimensions of null.stat matrix
# let's use smaller values than defaults so tests are faster
B <- 2 # default is: round( m * 10 / s )
# include a covariate for tests (a vector for individuals)
# draw continuous (in LFA, drawing bernoulli causes fitting problems I don't particularly want to deal with here)
# only perform these tests when qvalue is available (it isn't on some CRAN servers!)
if ( requireNamespace( "qvalue", quietly = TRUE ) ) {
# to prevent some qvalue-specific errors, ensure one p-value is 1
expect_error( pip( pvalue = pvalue, group = group[-1], verbose = FALSE ) )
prob <- pip( pvalue = pvalue, group = group, pi0 = pi0, verbose = FALSE )
# expect_true( all( prob <= 0 ) ) # not true!
prob <- pip( pvalue = pvalue, pi0 = pi0, verbose = FALSE )
# a rare error, ultimately in `qvalue::pi0est`, happens in the next line (due to small groups, not enough p-values for good pi0 estimates):
## Error (test-jackstraw.R:67:5): pip works
## Error: missing or infinite values in inputs are not allowed
##   9. jackstraw::pip(pvalue = pvalue, group = group, verbose = FALSE)
##  10. qvalue::lfdr(pvalue[group == i], ...) /home/viiia/docs/ochoalab/jackstraw/R/pip.R:45:16
##  12. stats::smooth.spline(lambda, pi0, df = smooth.df)
# so happens only when there's groups and pi0 has to be estimated
prob <- pip( pvalue = pvalue, group = group, verbose = FALSE )
prob <- pip( pvalue = pvalue, verbose = FALSE )
# (it could be an intercept only, but in practice we always pass objects constructed via model.matrix, so they're always matrices)
expect_error( RSS( dat = X, mod = 1:n ) )
rss_lm[ i ] <- sum( lm( X[ i, ] ~ LF1 )$residuals^2 )
# (not used by dependencies, but in theory supported)
# LF1 doesn't work as-is because the intercept gets added twice!
expect_error( FSTAT( dat = 1:10, LV = LV ) )
# `LV` has non-matching numbers of rows with `dat`
expect_error( FSTAT( dat = Xc, LV = LV[ -1, ] ) )
# since parammetric = FALSE, this only returns one element in its list
fstat_lm[ i ] <- summary( lm( X[ i, ] ~ LF1 ) )$fstatistic[1]
obj <- FSTAT( dat = Xc, LV = LV, covariate = covariate )
# successful run with ALV (should be identical to covariate, but some use cases specify both separately)
obj2 <- FSTAT( dat = Xc, LV = LV, ALV = covariate )
# NOTE: `parametric = TRUE` is not used in this package for any public code, so it's also not tested
#FSTAT(dat, LV, ALV = NULL, covariate = NULL, parametric = FALSE)
expect_error( permutationPA( dat = 1:10 ) )
obj <- permutationPA( dat = X, B = B, verbose = FALSE )
expect_equal( names( obj ), c('r', 'p') )
# a toy example of extremely structured data, where everybody is homozygous for one or the other allele
x <- c( rep.int( 0L, n2 ), rep.int( 2L, n2 ) )
# and there are no missing values in this case
# confirm that everything is in desired range
expect_true( all( y %in% c(0L, 1L, 2L) ) )
# for this example, the case without heterozygotes is so rare (0.5^10 = 9.8e-4) let's demand that there be at least one heterozygote
# and sum of genotypes equals `n` (in this case true by construction, not true in general)
# now try on the simulated random data (all rows)
# vectorize so number of tests isn't too extreme
Y <- t( apply( X, 1L, permute_alleles_from_geno ) )
# perform general tests, still allowing no missingness (not present in simulated data)
expect_true( all( Y %in% c(0L, 1L, 2L) ) )
# returning the same data is practically impossible if permutation is correct
# here test that the sums of input and output for each row match
expect_equal( rowSums( Y ), rowSums( X ) )
# we haven't performed tests with missingness in general, and the overall runtime is so high I don't think we want to do it broadly, but here let's just sprinkle random missingness (here there's absolute tolerance for bad things happening, like fixed rows, even rows full of NAs!)
X_miss[ sample( n*m, n*m*p_miss ) ] <- NA
Y_miss <- t( apply( X_miss, 1L, permute_alleles_from_geno ) )
# perform general tests, now allowing missingness!
expect_true( all( Y %in% c(0L, 1L, 2L), na.rm = TRUE ) )
expect_true( !all( Y == X, na.rm = TRUE ) )
expect_equal( rowSums( Y, na.rm = TRUE ), rowSums( X, na.rm = TRUE ) )
test_that( "empPvals handles NAs correctly", {
# `qvalue::empPvals` actually it doesn't, but it's easy to fix with a minor hack, wrapped around internal `empPvals`
expect_equal( length( pvals ), length( obs ) )
expect_true( all( is.na( pvals[ is.na( obs ) ] ) ) )
expect_true( !anyNA( pvals[ !is.na( obs ) ] ) )
test_jackstraw_return_val <- function ( obj, s, B, kmeans = FALSE ) {
# all jackstraw variants return basically the same thing
# kmeans object returns different names and order, unfortunately (otherwise equivalent though)
names_obj <- c('call', name_o, name_n, name_p)
names_obj <- c('call', name_p, name_o, name_n)
expect_true( !anyNA( p ) ) # in theory there can be NAs, they just don't arise in my simple examples
expect_true( all( p >= 0, na.rm = TRUE ) )
expect_true( all( p <= 1, na.rm = TRUE ) )
# NOTE: lengths of elements vary, not a useful test
expect_equal( length( null.stat ), s * B )
FUN <- function( x ) svd( x )$v[ , 1:d, drop = FALSE ]
# cause errors due to missing required data
# must provide all of dat = Xc, r = d, and FUN for a minimal successful run
expect_error( jackstraw_subspace( dat = Xc ) )
expect_error( jackstraw_subspace( r = d ) )
expect_error( jackstraw_subspace( FUN = FUN ) )
expect_error( jackstraw_subspace( dat = Xc, r = d ) )
expect_error( jackstraw_subspace( dat = Xc, FUN = FUN ) )
expect_error( jackstraw_subspace( r = d, FUN = FUN ) )
expect_error( jackstraw_subspace( dat = 1:10, r = d, FUN = FUN ) )
expect_error( jackstraw_subspace( dat = Xc, r = d, FUN = FUN, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_subspace( dat = Xc, r = d, FUN = FUN, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_subspace( dat = Xc, r = d, FUN = FUN, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_subspace( dat = Xc, r = d, FUN = FUN, covariate = rbind( covariate ) ) ) # transposed matrix
# make it silent so we can focus on problem messages
obj <- jackstraw_subspace( Xc, r = d, FUN = FUN, s = s, B = B, verbose = FALSE )
obj <- jackstraw_subspace( Xc, r = d, FUN = FUN, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_subspace( Xc, r = d, FUN = FUN, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_subspace( Xc, r = d, FUN = FUN, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_subspace( Xc, r = d, FUN = FUN, s = s, B = B, covariate = covariate, verbose = FALSE )
# must provide dat = Xc for a minimal successful run
expect_error( jackstraw_pca( dat = 1:10, r = d ) )
expect_error( jackstraw_pca( dat = Xc, r = d, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_pca( dat = Xc, r = d, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_pca( dat = Xc, r = d, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_pca( dat = Xc, r = d, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_pca( Xc, r = d, s = s, B = B, verbose = FALSE )
obj <- jackstraw_pca( Xc, r = d, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_pca( Xc, r = d, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_pca( Xc, r = d, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_pca( Xc, r = d, s = s, B = B, covariate = covariate, verbose = FALSE )
expect_error( jackstraw_rpca( dat = 1:10, r = d ) )
expect_error( jackstraw_rpca( dat = Xc, r = d, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_rpca( dat = Xc, r = d, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_rpca( dat = Xc, r = d, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_rpca( dat = Xc, r = d, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_rpca( Xc, r = d, s = s, B = B, verbose = FALSE )
obj <- jackstraw_rpca( Xc, r = d, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_rpca( Xc, r = d, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_rpca( Xc, r = d, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_rpca( Xc, r = d, s = s, B = B, covariate = covariate, verbose = FALSE )
expect_error( jackstraw_irlba( r = d ) )
expect_error( jackstraw_irlba( dat = 1:10, r = d ) )
expect_error( jackstraw_irlba( dat = Xc, r = d, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_irlba( dat = Xc, r = d, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_irlba( dat = Xc, r = d, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_irlba( dat = Xc, r = d, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_irlba( Xc, r = d, s = s, B = B, verbose = FALSE )
obj <- jackstraw_irlba( Xc, r = d, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_irlba( Xc, r = d, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_irlba( Xc, r = d, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_irlba( Xc, r = d, s = s, B = B, covariate = covariate, verbose = FALSE )
# must provide both dat = X and kmeans.dat for a minimal successful run
expect_error( jackstraw_kmeans( dat = Xc ) )
expect_error( jackstraw_kmeans( kmeans.dat = kmeans.dat ) )
expect_error( jackstraw_kmeans( dat = 1:10, kmeans.dat = kmeans.dat ) )
expect_error( jackstraw_kmeans( dat = Xc, kmeans.dat = kmeans.dat, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_kmeans( dat = Xc, kmeans.dat = kmeans.dat, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_kmeans( dat = Xc, kmeans.dat = kmeans.dat, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_kmeans( dat = Xc, kmeans.dat = kmeans.dat, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_kmeans( Xc, kmeans.dat = kmeans.dat, s = s, B = B, verbose = FALSE )
test_jackstraw_return_val( obj, s, B, kmeans = TRUE )
obj <- jackstraw_kmeans( Xc, kmeans.dat = kmeans.dat, s = 1, B = B, verbose = FALSE )
test_jackstraw_return_val( obj, 1, B, kmeans = TRUE )
obj <- jackstraw_kmeans( Xc, kmeans.dat = kmeans.dat, s = s, B = 1, verbose = FALSE )
test_jackstraw_return_val( obj, s, 1, kmeans = TRUE )
obj <- jackstraw_kmeans( Xc, kmeans.dat = kmeans.dat, s = 1, B = 1, verbose = FALSE )
test_jackstraw_return_val( obj, 1, 1, kmeans = TRUE )
obj <- jackstraw_kmeans( Xc, kmeans.dat = kmeans.dat, s = s, B = B, covariate = covariate, verbose = FALSE )
kmeans.dat <- ClusterR::KMeans_rcpp( Xc, clusters = d )
expect_error( jackstraw_kmeanspp( dat = Xc ) )
expect_error( jackstraw_kmeanspp( kmeans.dat = kmeans.dat ) )
expect_error( jackstraw_kmeanspp( dat = 1:10, kmeans.dat = kmeans.dat ) )
expect_error( jackstraw_kmeanspp( dat = Xc, kmeans.dat = kmeans.dat, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_kmeanspp( dat = Xc, kmeans.dat = kmeans.dat, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_kmeanspp( dat = Xc, kmeans.dat = kmeans.dat, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_kmeanspp( dat = Xc, kmeans.dat = kmeans.dat, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_kmeanspp( Xc, kmeans.dat = kmeans.dat, s = s, B = B, verbose = FALSE )
obj <- jackstraw_kmeanspp( Xc, kmeans.dat = kmeans.dat, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_kmeanspp( Xc, kmeans.dat = kmeans.dat, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_kmeanspp( Xc, kmeans.dat = kmeans.dat, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_kmeanspp( Xc, kmeans.dat = kmeans.dat, s = s, B = B, covariate = covariate, verbose = FALSE )
MiniBatchKmeans.output <- ClusterR::MiniBatchKmeans( Xc, clusters = d, batch_size = batch_size )
# must provide both dat = X and MiniBatchKmeans.output for a minimal successful run
expect_error( jackstraw_MiniBatchKmeans( dat = Xc ) )
expect_error( jackstraw_MiniBatchKmeans( MiniBatchKmeans.output = MiniBatchKmeans.output ) )
expect_error( jackstraw_MiniBatchKmeans( dat = 1:10, MiniBatchKmeans.output = MiniBatchKmeans.output ) )
expect_error( jackstraw_MiniBatchKmeans( dat = Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_MiniBatchKmeans( dat = Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_MiniBatchKmeans( dat = Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_MiniBatchKmeans( dat = Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_MiniBatchKmeans( Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, batch_size = batch_size, s = s, B = B, verbose = FALSE )
obj <- jackstraw_MiniBatchKmeans( Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, batch_size = batch_size, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_MiniBatchKmeans( Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, batch_size = batch_size, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_MiniBatchKmeans( Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, batch_size = batch_size, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_MiniBatchKmeans( Xc, MiniBatchKmeans.output = MiniBatchKmeans.output, batch_size = batch_size, s = s, B = B, covariate = covariate, verbose = FALSE )
# must provide both dat = X and pam.dat for a minimal successful run
expect_error( jackstraw_pam( dat = Xc ) )
expect_error( jackstraw_pam( pam.dat = pam.dat ) )
expect_error( jackstraw_pam( dat = 1:10, pam.dat = pam.dat ) )
expect_error( jackstraw_pam( dat = Xc, pam.dat = pam.dat, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_pam( dat = Xc, pam.dat = pam.dat, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_pam( dat = Xc, pam.dat = pam.dat, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_pam( dat = Xc, pam.dat = pam.dat, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_pam( Xc, pam.dat = pam.dat, s = s, B = B, verbose = FALSE )
# NOTE: removed s=1 cases because, though they are fine often, sometimes there are sigularity errors that are not worth debugging for such small toy cases tested here
##     ###ERROR: Error in `solve.default(t(mod) %*% mod)`: Lapack routine dgesv: system is exactly singular: U[2,2] = 0
##     obj <- jackstraw_pam( Xc, pam.dat = pam.dat, s = 1, B = B, verbose = FALSE )
## test_jackstraw_return_val( obj, 1, B, kmeans = TRUE )
obj <- jackstraw_pam( Xc, pam.dat = pam.dat, s = s, B = 1, verbose = FALSE )
##     obj <- jackstraw_pam( Xc, pam.dat = pam.dat, s = 1, B = 1, verbose = FALSE )
## test_jackstraw_return_val( obj, 1, 1, kmeans = TRUE )
obj <- jackstraw_pam( Xc, pam.dat = pam.dat, s = s, B = B, covariate = covariate, verbose = FALSE )
# 4 arguments are required for a successful run
expect_error( jackstraw_cluster( dat = Xc ) )
expect_error( jackstraw_cluster( k = d ) )
expect_error( jackstraw_cluster( cluster = cluster ) )
expect_error( jackstraw_cluster( centers = centers ) )
expect_error( jackstraw_cluster( dat = Xc, k = d ) )
expect_error( jackstraw_cluster( dat = Xc, cluster = cluster ) )
expect_error( jackstraw_cluster( dat = Xc, centers = centers ) )
expect_error( jackstraw_cluster( k = d, cluster = cluster ) )
expect_error( jackstraw_cluster( k = d, centers = centers ) )
expect_error( jackstraw_cluster( cluster = cluster, centers = centers ) )
expect_error( jackstraw_cluster( k = d, cluster = cluster, centers = centers ) )
expect_error( jackstraw_cluster( dat = Xc, cluster = cluster, centers = centers ) )
expect_error( jackstraw_cluster( dat = Xc, k = d, centers = centers ) )
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster ) )
expect_error( jackstraw_cluster( dat = 1:10, k = d, cluster = cluster, centers = centers ) )
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster[-1], centers = centers ) )
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers[ -1, ] ) )
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers[ , -1 ] ) )
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, algorithm = FUN, s = s, B = B, verbose = FALSE )
obj <- jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, algorithm = FUN, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, algorithm = FUN, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, algorithm = FUN, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_cluster( dat = Xc, k = d, cluster = cluster, centers = centers, algorithm = FUN, s = s, B = B, covariate = covariate, verbose = FALSE )
file_tmp <- tempfile( 'test-jackstraw', fileext = '.bed' )
genio::write_bed( file_tmp, X, verbose = FALSE )
# read it this way, specifying dimensions, because there's no BIM/FAM files
X_BM <- BEDMatrix::BEDMatrix( file_tmp, n = n, p = m )
# NOTE: BEDMatrix and genio are both dependencies
# recall s==2, let's leave it that way for now
obj <- jackstraw_BEDMatrix( dat = X_BM, random_s, m_chunk = m_chunk )
expect_equal( names( obj ), c('dat_full', 'dat_rand', 'file_full', 'file_rand') )
expect_true( class( obj$dat_full ) == 'BEDMatrix' )
expect_true( class( obj$dat_rand ) == 'BEDMatrix' )
expect_equal( class( obj$file_full ), 'character' )
expect_equal( class( obj$file_rand ), 'character' )
# the non-random loci should match the original data
# (NOTE: BEDMatrix data has to be transposed and dimnames removed)
X_orig_obs <- t( obj$dat_full[ , -random_s ] )
# the random data in the full matrix should equal the random matrix
# here they're both BEDMatrix so no transposition is required
# however, `dat_rand` lists rows in original order, which may disagree with `random_s` (random order), so let's sort here to get matching orders
X_rand_obs <- obj$dat_full[ , random_s_sorted, drop = FALSE ]
# need sorted version for some tests (not needed if s=1 but meh)
# the following tests require the Bioconductor `lfa` package only
if ( requireNamespace( "lfa", quietly = TRUE ) ) {
expect_error( efron_Rsq_snp( snp = xi ) )
expect_error( efron_Rsq_snp( p1 = pi ) )
r2 <- efron_Rsq_snp( snp = xi, p1 = pi )
# the basics of what this R^2 should be like
expect_error( mcfadden_Rsq_snp( snp = xi ) )
expect_error( mcfadden_Rsq_snp( p1 = p1 ) )
expect_error( mcfadden_Rsq_snp( p0 = p0 ) )
expect_error( mcfadden_Rsq_snp( snp = xi, p1 = p1 ) )
expect_error( mcfadden_Rsq_snp( snp = xi, p0 = p0 ) )
expect_error( mcfadden_Rsq_snp( p1 = p1, p0 = p0 ) )
r2 <- mcfadden_Rsq_snp( snp = xi, p1 = p1, p0 = p0 )
expect_error( pseudo_Rsq( LF_alt = LF1 ) )
expect_error( pseudo_Rsq( as.vector( X ), LF1 ) )
expect_error( pseudo_Rsq( X, as.vector( LF1 ) ) )
# LF_null is set to default (intercept only)
expect_error( efron_Rsq( as.vector( X ), LF1 ) )
expect_error( efron_Rsq( X, as.vector( LF1 ) ) )
# these functions require the Bioconductor package `gcatest`
if ( requireNamespace( "gcatest", quietly = TRUE ) ) {
# the following tests require the Bioconductor `lfa` package
# define the function to pass to `jackstraw_lfa`!  Uses global `d`
# must provide all of (dat = X, r = d, FUN = FUN) for a minimal successful run
expect_error( jackstraw_lfa( dat = X ) )
expect_error( jackstraw_lfa( FUN = FUN ) )
expect_error( jackstraw_lfa( dat = X, r = d ) )
expect_error( jackstraw_lfa( dat = X, FUN = FUN ) )
expect_error( jackstraw_lfa( r = d, FUN = FUN ) )
expect_error( jackstraw_lfa( dat = 1:10, r = d, FUN = FUN ) )
expect_error( jackstraw_lfa( dat = X, r = d, FUN = FUN, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_lfa( dat = X, r = d, FUN = FUN, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_lfa( dat = X, r = d, FUN = FUN, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_lfa( dat = X, r = d, FUN = FUN, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_lfa( X, r = d, FUN = FUN, s = s, B = B, verbose = FALSE )
obj <- jackstraw_lfa( X, r = d, FUN = FUN, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_lfa( X, r = d, FUN = FUN, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_lfa( X, r = d, FUN = FUN, s = 1, B = 1, verbose = FALSE )
# this comparison succeeds most of the time, but not 100% of the time
# bad fits cause errors randomly, which are rare but over 100 loci it gets less rare
obj <- jackstraw_lfa( X, r = d, FUN = FUN, s = s, B = B, covariate = covariate, verbose = FALSE )
# test version without default allele-level permutation!
obj <- jackstraw_lfa( X, r = d, FUN = FUN, s = s, B = B, permute_alleles = FALSE, verbose = FALSE )
test_that( "jackstraw_lfa works with BEDMatrix", {
obj <- jackstraw_lfa( X_BM, r = d, FUN = FUN, s = s, B = B, verbose = FALSE )
obj <- jackstraw_lfa( X_BM, r = d, FUN = FUN, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_lfa( X_BM, r = d, FUN = FUN, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_lfa( X_BM, r = d, FUN = FUN, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_lfa( X_BM, r = d, FUN = FUN, s = s, B = B, covariate = covariate, verbose = FALSE )
# run alstructure tests only if optional package is available
# define function to pass (uses global `d`)
FUN <- function(x) t( alstructure(x, d_hat = d)$Q_hat )
expect_error( jackstraw_alstructure( dat = X ) )
expect_error( jackstraw_alstructure( r = d ) )
expect_error( jackstraw_alstructure( FUN = FUN ) )
expect_error( jackstraw_alstructure( dat = X, r = d ) )
expect_error( jackstraw_alstructure( dat = X, FUN = FUN ) )
expect_error( jackstraw_alstructure( r = d, FUN = FUN ) )
expect_error( jackstraw_alstructure( dat = 1:10, r = d, FUN = FUN ) )
expect_error( jackstraw_alstructure( dat = X, r = d, FUN = FUN, covariate = 1 ) ) # scalar is bad
expect_error( jackstraw_alstructure( dat = X, r = d, FUN = FUN, covariate = covariate[-1] ) ) # length off by 1, vector
expect_error( jackstraw_alstructure( dat = X, r = d, FUN = FUN, covariate = cbind( covariate[-1] ) ) ) # length off by 1, matrix
expect_error( jackstraw_alstructure( dat = X, r = d, FUN = FUN, covariate = rbind( covariate ) ) ) # transposed matrix
obj <- jackstraw_alstructure( X, r = d, FUN = FUN, s = s, B = B, verbose = FALSE )
obj <- jackstraw_alstructure( X, r = d, FUN = FUN, s = 1, B = B, verbose = FALSE )
obj <- jackstraw_alstructure( X, r = d, FUN = FUN, s = s, B = 1, verbose = FALSE )
obj <- jackstraw_alstructure( X, r = d, FUN = FUN, s = 1, B = 1, verbose = FALSE )
obj <- jackstraw_alstructure( X, r = d, FUN = FUN, s = s, B = B, covariate = covariate, verbose = FALSE )
invisible( suppressWarnings( file.remove( file_tmp ) ) )

# this has to be at least this large or `pip` (really `qvalue::pi0est`) fails in stupid ways
# in theory there can be NAs, they just don't arise in my simple examples