Adjust an SCV value for the bias arising when it is calculated from unbiased estimates of mean and variance.

Description

Assume that a small sample of i.i.d. random variables from a negative binomial distribution is given, and you have obtained unbiased estimates of mean and raw variance. Then, a new bias is introduced when the squared coefficient of variation (SCV, a.k.a. dispersion) is calculated from these unbiased estimates by dividing the raw variance by the square of the mean. This bias can be calculated by numerical simulation and a pre-calculated adjustment table (or rather a fit through tabulated values) is supplied with the package. The present function uses this to remove the bias from a raw SCV estimate.

This function is used internally in nbinomTest. You will rarely need to call it directly.

Usage

adjustScvForBias(scv, nsamples)

Arguments

Value

an unbiased estimate of the raw SCV

Author(s)

Simon Anders

Examples

   true_mean <- 100
   true_scv <- .1
   nsamples <- 3
   res <- replicate( 1000, {
      mySample <- rnbinom( nsamples, mu=true_mean, size=1/true_scv )
      mu_est <- mean( mySample )
      raw_var_est <- var( mySample ) - mean( mySample )
      raw_scv_est <- raw_var_est / mu_est^2
      unbiased_raw_scv_est <- adjustScvForBias( raw_scv_est, 4 )
      c( raw_scv_est = raw_scv_est, unbiased_raw_scv_est = unbiased_raw_scv_est ) } )
   rowMeans( res )

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(DESeq)
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: locfit
locfit 1.5-9.1 	 2013-03-22
Loading required package: lattice
    Welcome to 'DESeq'. For improved performance, usability and
    functionality, please consider migrating to 'DESeq2'.
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/DESeq/adjustScvForBias.Rd_%03d_medium.png", width=480, height=480)
> ### Name: adjustScvForBias
> ### Title: Adjust an SCV value for the bias arising when it is calculated
> ###   from unbiased estimates of mean and variance.
> ### Aliases: adjustScvForBias
> 
> ### ** Examples
> 
>    true_mean <- 100
>    true_scv <- .1
>    nsamples <- 3
>    res <- replicate( 1000, {
+       mySample <- rnbinom( nsamples, mu=true_mean, size=1/true_scv )
+       mu_est <- mean( mySample )
+       raw_var_est <- var( mySample ) - mean( mySample )
+       raw_scv_est <- raw_var_est / mu_est^2
+       unbiased_raw_scv_est <- adjustScvForBias( raw_scv_est, 4 )
+       c( raw_scv_est = raw_scv_est, unbiased_raw_scv_est = unbiased_raw_scv_est ) } )
>    rowMeans( res )
         raw_scv_est unbiased_raw_scv_est 
           0.1025996            0.1087198 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>