genomic.autocorr (version 1.0-1)

block.glm: block.glm

Description

Regression models for genomic data often assume there is independence between neighbouring genomic elements when, in reality, there is spatial dependence. This function implements a block bootstrap method for estimating correct variances of parameter estimates.

Usage

block.glm(f.lhs, f.rhs, data, order.by = NULL, strat.by = NULL,
  block.size = 20, B = 200, ...)

Arguments

f.lhs

character vector, left hand side of a formula, the model(s) to be fit will be defined by `paste(f.lhs, f.rhs, sep=" ~ ")`

f.rhs

character string, right hand side of a formula

data

data.table containing the columns referred to in f.lhs and f.rhs

order.by

if not `NULL`, the name of a column in `data` on which it should be sorted

strat.by

if not `NULL`, the name of a column in `data` on which it should be stratified before block sampling. Eg, if you are considering genomic data, you should stratify by chromosome as there should be no spatial correlation between chromosomes

block.size

size of blocks of contiguous observations that will be sampled for bootstrap estimation of variance of parameter estimates

B

number of bootstrap estimates

...

other arguments passed to `glm()` (eg `family="binomial"`)

Value

data.table giving the estimated effect ("beta") of each item in f.rhs on each item in f.lhs, together with block bootstrap estimates of confidence interval (beta.025, beta.975) and standard error (se.beta) and the number of bootstraps on which those estimates are based.

Details

Note that this function uses `mclapply` to parallelise the bootstrapping. Please set `mc.cores` to something sensible, eg options(mc.cores=10) if you have 10 cores.

Examples

# NOT RUN {
## simulate data with 10 repeated observations in a row - ie there
## should be autocorrelation only within windows <= 10
library(data.table)
data <- genomic.autocorr:::.sim.data(beta=0.2) 

## suppose we ignored the autocorrelation and look at the
## confidence interval for the effect of x on y1
r1<-summary(glm(y1 ~ x, data=data))$coefficients
r1

## if we know the block structure, as here, we can see the
## confidence interval is (inappropriately) much tighter than
## if we used just independent observations
r2<-summary(glm(y1 ~ x, data=data[!duplicated(name),]))$coefficients
r2

## use block bootstrap - x should only have a significant effect
## on y1 and the confidence interval around its effect should be
## closer to r2, above
r <- block.glm(f.lhs=c("y0","y1"), f.rhs="x",data=data,block.size=20,B=200)
r

## compare the block bootstrap and model based confidence intervals for x on y1
results <- rbind(c(r1[2,1], r1[2,1]-1.96*r1[2,2], r1[2,1]+1.96*r1[2,2]),
c(r2[2,1], r2[2,1]-1.96*r2[2,2], r2[2,1]+1.96*r2[2,2]),
as.numeric(r[4,.(beta,beta.025,beta.975)]))
dimnames(results) <- list(c("standard, ignore blocked","standard, independent obs","bootstrap"),
c("beta","LCI","UCI"))
results

with(as.data.frame(results), {
plot(1:nrow(results), beta,ylim=c(min(c(-0.01,LCI)),max(UCI)),axes=FALSE,xlab="Method",
main="Comparison of confidence intervals around coefficient estimates")
segments(x0=1:nrow(results),y0=LCI,y1=UCI)
abline(h=c(0,0.2),lty="dotted")
axis(1,1:nrow(results),rownames(results))
axis(2)
text(x=c(3,3),y=c(0,0.2),labels=c("null","true"),adj=c(1.1,0))
box()
})

# }