MahShaaban / pcr

Quality assessing, analyzing and testing the statistical significance of real-time quantitative PCR data
https://CRAN.R-project.org/package=pcr
GNU General Public License v3.0
27 stars 8 forks source link

made tidy df optional, default #9

Closed MahShaaban closed 5 years ago

MahShaaban commented 5 years ago

The current return value of pcr_test is a tidy data.frame. I added an argument tidy with default TRUE so it returns the same. When tidy = FALSE the return value, is a list of the original test object for each gene; lm when the required test = 'lm' and htest when test = 't.test' or test = 'wilcox.test'. This enables all sort of things to be done to the original test object. Below is an example with calling anova(), residuals() and fitted on the lm() object.

library(pcr)

# load ct values
fl <- system.file('extdata', 'ct4.csv', package = 'pcr')
ct4 <- readr::read_csv(fl)
#> Parsed with column specification:
#> cols(
#>   ref = col_double(),
#>   target = col_double()
#> )

# make group variable
group <- rep(c('control', 'treatment'), each = 12)

# test using pcr_test and return a list of lm object
res <- pcr_test(ct4,
                group_var = group,
                reference_gene = 'ref',
                reference_group = 'control',
                test = 'lm',
                tidy = FALSE)

# call anova on the lm object
anova(res$target)
#> Analysis of Variance Table
#> 
#> Response: x
#>           Df Sum Sq Mean Sq F value    Pr(>F)    
#> group_var  1 2.8139 2.81391  27.643 2.829e-05 ***
#> Residuals 22 2.2395 0.10179                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# extract residuals
residuals(res$target)
#>            1            2            3            4            5 
#> -0.271708333 -0.234108333 -0.133808333  0.955891667  0.129991667 
#>            6            7            8            9           10 
#> -0.149808333 -0.338808333  0.116791667 -0.279708333 -0.095108333 
#>           11           12           13           14           15 
#>  0.005691667  0.294691667  0.378916667 -0.021883333  0.379316667 
#>           16           17           18           19           20 
#> -0.415883333 -0.294083333 -0.078883333  0.191616667  0.340916667 
#>           21           22           23           24 
#> -0.234983333  0.110616667 -0.174183333 -0.181483333

# extract fitted values
fitted(res$target)
#>        1        2        3        4        5        6        7        8 
#> 3.640408 3.640408 3.640408 3.640408 3.640408 3.640408 3.640408 3.640408 
#>        9       10       11       12       13       14       15       16 
#> 3.640408 3.640408 3.640408 3.640408 2.955583 2.955583 2.955583 2.955583 
#>       17       18       19       20       21       22       23       24 
#> 2.955583 2.955583 2.955583 2.955583 2.955583 2.955583 2.955583 2.955583

# make diagnostic plots
## qq-plot of the residuals
qqnorm(residuals(res$target))


## fitted vs residuals
plot(x = fitted(res$target),
     y = residuals(res$target))