/
test_degree.R
77 lines (72 loc) · 2.82 KB
/
test_degree.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#' test_degree
#'
#' Test if a degree of relatedness is consistent with an observed PMR
#'
#' @param in_tibble a tibble that is the output of the callRelatedness() function.
#' @param row either the row number or pair name for which the posterior distribution is to be plotted.
#' @param degree the degree of relatedness to be tested.
#' @param verbose a logical (boolean) for whether all test output should be printed to screen.
#'
#' @return the associated p-value for the test
#' @export
#'
#' @examples
#' test_degree(relatedness_example, 1, 1)
test_degree <- function(in_tibble,row,degree,verbose=TRUE){
# Test that the in_tibble is of the correct form
if(nrow(in_tibble)==0){
stop('The input tibble is empty')
}
col_names_check <- c('row','pair','relationship','pmr','sd','mismatch','nsnps','ave_rel','Same_Twins','First_Degree','Second_Degree','Unrelated')
if((ncol(in_tibble)!=12)|(!all(names(in_tibble)==col_names_check))){
stop(paste0('Input tibble/data.frame must have 4 columns named: ',paste0(col_names_check,collapse=', '),'.'))
}
# check row is a numeric or char variable, and if char, find the row.
if(is.character(row)){
rowtmp <- row
row <- which(row==in_tibble$pair)
if(length(row)==0){
stop(paste0('There is no pair called: ',rowtmp))
}
}else if(is.numeric(row)&(row>nrow(in_tibble))){
stop('row number requested was higher than the number of rows in input tibble!')
}else if(!is.numeric(row)){
stop('row must be a number or a pair name.')
}
if(!is.numeric(degree)){
stop('Degree must be numeric.')
}
# if(degree>10|degree<0|(degree%%1!=0)){ # Removed constraint that integer
if(degree>10|degree<0){
stop('Degree must be a value k such that 0<k<10.')
}
M <- in_tibble$ave_rel[row]
pk <- (1-0.5^(degree+1))*M
N <- in_tibble$nsnps[row]
x <- in_tibble$mismatch[row]
pair <- in_tibble$pair[row]
# pval <- 2*stats::pbinom(x,N,pk,ifelse(x>pk*N,F,T))
pval <- min(
c(
2*stats::pbinom(x,N,pk,T),
2*stats::pbinom(x,N,pk,F)
)
)
# CI <- x/N+c(-1,1)*1.96*sqrt(x/N*(1-x/N)/N)
pvaltext <- ifelse(pval<0.0001,
sprintf('p-value : %s\n',format(pval,scientific=T)),
sprintf('p-value : %.4f\n',pval))
degree_string <- dplyr::case_when(degree==1 ~ "1st",
degree==2 ~ "2nd",
degree==3 ~ "3rd",
T ~ paste0(degree,"th"))
if(verbose){
cat(sprintf('Testing H0 : "%s" are %s-degree relatives.\n',pair,degree_string))
cat(sprintf('Expected PMR : %.4f\n',pk))
cat(sprintf('Observed PMR : %.4f\n',x/N))
cat(sprintf('Estimated degree : %.4f\n',convertP(x/N,M)))
cat(pvaltext)
cat(ifelse(pval>=0.05,'Decision : Retain H0\n','Decision : Reject H0\n'))
}
return(pval)
}