summary.bcp {bcp} | R Documentation |
summary method for class ``bcp''.
summary.bcp(object, quantiles=c(0.025, 0.25, 0.5, 0.75, 0.975), digits = max(3, .Options$digits - 3), ...)
object |
the result of a call to `bcp'. |
quantiles |
an optional vector of quantiles. |
digits |
the number of digits displayed in the summary statistics. |
... |
additional arguments. |
The code for `summary.bcp' was taken from the `summary.mcmc' function found in the coda package (Plummer et al., 2006). The function returns the posterior probability of a change point for each position, the posterior means and their quantiles over the mcmc iterations, along with the standard deviation and standard error of the mean.
Chandra Erdman
Bai J., Perron P. (2003), Computation and Analysis of Multiple Structural Change Models, Journal of Applied Econometrics, 18, 1-22. url: http: //qed.econ.queensu.ca/jae/2003-v18.1/bai-perron/.
Daniel Barry and J. A. Hartigan (1993), A Bayesian Analysis for Change Point Problems, Journal of The American Statistical Association, 88, 309-19.
Olshen, A. B., Venkatraman, E. S., Lucito, R., Wigler, M. (2004), Circular binary segmentation for the analysis of array-based DNA copy number data, Biostatistics, 5, 557-572. url: http://www.bioconductor.org/repository/release1.5/package/html/DNAcopy.html.
Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines (2006), The coda Package, version 0.10-7, CRAN: The Comprehensive R Network.
Snijders et al. (2001), Assembly of microarrays for genome-wide measurement of DNA copy number, Nature Genetics, 29, 263-264.
Achim Zeileis, Friedrich Leisch, Bruce Hansen, Kurt Hornik, Christian Kleiber (2006), The strucchange Package, version 1.3-1, CRAN: The Comprehensive R Network.
`bcp', `print.bcp', and `plot.bcp'.
##### A random sample from a few normal distributions ##### testdata <- c(rnorm(20), rnorm(10, 5, 1), rnorm(20)) bcp.0 <- bcp(testdata) summary.bcp(bcp.0) plot.bcp(bcp.0) ##### Coriell chromosome 11 ##### data(coriell) chrom11 <- na.omit(coriell$Coriell.05296[coriell$Chromosome==11]) n <- length(chrom11) bcp.11 <- bcp(chrom11[1:n]) summary.bcp(bcp.11) plot.bcp(bcp.11) # to see bcp and Circular Binary Segmentation results run: if(require("DNAcopy")) { bcp.11$rhos[,ncol(bcp.11$rhos)] <- 0 changepoint.freq <- apply(bcp.11$rhos[bcp.11$burnin:(bcp.11$burnin+bcp.11$burnin),1:dim(bcp.11$results)[2]],2,mean) cbs <- segment(CNA(chrom11[1:n], rep(1, n), 1:n), verbose = 0) cbs.ests <- rep(unlist(cbs$output[6]), unlist(cbs$output[5])) op <- par(mfrow=c(2,1),col.lab="black",col.main="black") plot(1:n, bcp.11$posterior.mean, type="l", xlab="Location", ylab="Posterior Mean", main="Posterior Means") lines(cbs.ests, col="red") points(chrom11) plot(1:n, changepoint.freq, type="l", xlab="Location", ylab="Posterior Probability of a Change", main="Change Point Locations") for(i in 1:(dim(cbs$output)[1]-1)) abline(v=cbs$output$loc.end[i], col="red") par(op) } else { cat("DNAcopy is not loaded") } ##### RealInt ##### data("RealInt") bcp.ri <- bcp(as.vector(RealInt)) summary.bcp(bcp.ri) plot.bcp(bcp.ri) # to see bcp and Bai and Perron results run: if(require("strucchange")) { bcp.ri$rhos[,ncol(bcp.ri$rhos)] <- 0 changepoint.freq <- apply(bcp.ri$rhos[bcp.ri$burnin:(bcp.ri$burnin+bcp.ri$burnin),1:dim(bcp.ri$results)[2]],2,mean) bp <- breakpoints(RealInt ~ 1, h = 2)$breakpoints rho <- rep(0, length(RealInt)) rho[bp] <- 1 b.num<-1 + c(0,cumsum(rho[1:(length(rho)-1)])) bp.mean <- unlist(lapply(split(RealInt,b.num),mean)) bp.ri <- rep(0,length(RealInt)) for(i in 1:length(bp.ri)) bp.ri[i] <- bp.mean[b.num[i]] op <- par(mfrow=c(2,1),col.lab="black",col.main="black") xax <- seq(1961, 1987, length=103) plot(xax, bcp.ri$posterior.mean, type="l", xlab="Time", ylab="Mean", main="Posterior Means") lines(xax, bp.ri, col="blue") points(RealInt) plot(xax, changepoint.freq, type="l", xlab="Time", ylab="Posterior Probability", main="Posterior Probability of a Change") for(i in 1:length(bp.ri)) abline(v=xax[bp[i]], col="blue") par(op) } else { cat("strucchange is not loaded") }