Showing posts with label R. Show all posts
Showing posts with label R. Show all posts

Thursday, April 10, 2014

R Contour / Density Scatter plot

http://stats.stackexchange.com/questions/31726/scatterplot-with-contour-heat-overlay

n = 10000;
a = as.matrix( rnorm(n) );
dim(a) =  c( n/2, 2 );
z = kde2d( a[,1], a[,2], n=50 );
plot( a[,1], a[,2], cex=0.5 );
contour( add=T, z, col=1:11 );

Tuesday, March 18, 2014

Multivariate analysis

http://little-book-of-r-for-multivariate-analysis.readthedocs.org/en/latest/src/multivariateanalysis.html

Multivariate Analysis

This booklet tells you how to use the R statistical software to carry out some simple multivariate analyses, with a focus on principal components analysis (PCA) and linear discriminant analysis (LDA).

PCA How To
http://psych.colorado.edu/wiki/lib/exe/fetch.php?media=labs:learnr:emily_-_principal_components_analysis_in_r:pca_how_to.pdf 

So that’s it. To do PCA, all you have to do is follow these steps:
1. Get X in the proper form. This will probably mean subtracting off the means of
each row. If the variances are significantly different in your data, you may also
wish to scale each row by dividing by its standard deviation to give the rows a
uniform variance of 1 (the subtleties of how this affects your analysis are beyond
the scope of this paper, but in general, if you have significantly different scales in
your data it’s probably a good idea).
2. Calculate A=XXT
3. Find the eigenvectors of A and stack them to make the matrix P.
4. Your new data is PX, the new variables (a.k.a. principal components) are the rows
of P.
5. The variance for each principal component can be read off the diagonal of the
covariance matrix.

# Obtain data in a matrix
Xoriginal=t(as.matrix(recorded.data))
# Center the data so that the mean of each row is 0
rm=rowMeans(Xoriginal)
X=Xoriginal-matrix(rep(rm, dim(X)[2]), nrow=dim(X)[1])
# Calculate P
A=X %*% t(X)
E=eigen(A,TRUE)
P=t(E$vectors)
# Find the new data and standard deviations of the principal components
newdata = P %*% X
sdev = sqrt(diag((1/(dim(X)[2]-1)* P %*% A %*% t(P))))

Going back to the derivation of PCA, we have N=PX, where N is our new data. Since we
know that P-1=PT , it is easy to see that X=PT N. Thus, if we know P and N, we can easily
recover X. This is useful, because if we choose to throw away some of the smaller
components—which hopefully are just noise anyway—N is a smaller dataset than X. But
we can still reconstruct data that is almost the same as X.

Run prcomp() again, but this time include the option tol=0.1. What is returned will be
any principal components whose standard deviation is greater than 10% of the standard
deviation of the first principal component. In this case, the first two components are
returned.

pr=prcomp(recorded.data)
pr
plot(pr)
barplot(pr$sdev/pr$sdev[1])
pr2=prcomp(recorded.data, tol=.1)
plot.ts(pr2$x)
quartz(); plot.ts(intensities)
quartz(); plot.ts(recorded.data)
quartz(); plot.ts(cbind(-1*pr2$x[,1],pr2$x[,2]))

You can see that od, which is the reconstruction of X  from when no principal components were discarded, is identical to the recorded.data.  od2 is the reconstruction of X from only two principal components.

Because prcomp() works with variables in columns instead of rows as in the derivation
above, the required transformation is X=NPT
, or in R syntax, X=pr$x %*% t(pr$rotation).
Run the code at the top of the page

od=pr$x %*% t(pr$rotation)
od2=pr2$x %*% t(pr2$rotation)
quartz(); plot.ts(recorded.data)
quartz(); plot.ts(od)
quartz(); plot.ts(od2)

Thursday, January 23, 2014

Software carpentry

http://software-carpentry.org/v4/index.html

Who We Are

Our volunteers teach basic software skills to researchers in science, engineering, and medicine. Founded in 1998, we are now part of the Mozilla Science Lab.

What We Do

We run bootcamps all over the world, and provide open access material for self-paced instruction. We also run a training program for people who'd like to help us teach.

How To Help

Like all volunteer organizations, we depend on you to help us help others. You can host a bootcamp, help create new teaching materials, or improve the tools we use.

Tuesday, January 21, 2014

Taking R to the Limit (High Performance Computing in R)

http://www.slideshare.net/bytemining/r-hpc

bigmemory - it is ideal for problems involving the analysis in R for manageable subsets of the data, or when an analysis is conducted mostly in C++

the "big" family

biganalytics, synchronicity, bigtabulate, big.matrix, bigalgebra, bigvideo, shared.big.matrix, filebacked.big.matrix, bigsplit

linear models: biglm.big.matrix

mwhich

ff - "fast access files" - file-based access to datasets that cannot fit in memory

data.table

mapReduce - apply(map(data), reduce)

HadoopStreaming

Wednesday, January 15, 2014

data.table: Extension of data.frame for fast indexing, fast ordered joins, fast assignment, fast grouping and list columns

http://cran.r-project.org/web/packages/data.table/index.html

Enhanced data.frame. Fast indexing, fast ordered joins, fast assignment by reference, fast grouping and list columns in a short and flexible syntax. i and j may be expressions of column names directly, for faster development. Example: X[Y] is a fast join for large data.

http://datatable.r-forge.r-project.org/ 

Friday, October 11, 2013

Surrogate variable analysis

http://www.bioconductor.org/packages/2.12/bioc/vignettes/sva/inst/doc/sva.pdf

The sva package contains functions for removing batch e ects and other unwanted variation
in high-throughput experiments. Speci cally, the sva package contains functions for identifying and building surrogate variables for high-dimensional data sets. Surrogate variables
are covariates constructed directly from high-dimensional data (like gene expression/RNA
sequencing/methylation/brain imaging data) that can be used in subsequent analyses to
adjust for unknown, unmodeled, or latent sources of noise.

Thursday, August 29, 2013

Finding local extrema of a density function using splines

http://stats.stackexchange.com/questions/30750/finding-local-extrema-of-a-density-function-using-splines

require(graphics)
#some data
d <- density(faithful$eruptions, bw = "sj")

#make it a time series
ts_y<-ts(d$y)

#calculate turning points (extrema)
require(pastecs)
tp<-turnpoints(ts_y)
#plot
plot(d)
points(d$x[tp$tppos],d$y[tp$tppos],col="red")

Monday, June 17, 2013

pheatmap

http://rgm3.lab.nig.ac.jp/RGM/r_function?p=pheatmap&f=pheatmap

library(pheatmap)

# Generate some data
test = matrix(rnorm(200), 20, 10)
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
colnames(test) = paste("Test", 1:10, sep = "")
rownames(test) = paste("Gene", 1:20, sep = "")


# Generate column annotations
annotation = data.frame(Var1 = factor(1:10 %% 2 == 0, labels = c("Class1", "Class2")), Var2 = 1:10)
annotation$Var1 = factor(annotation$Var1, levels = c("Class1", "Class2", "Class3"))
rownames(annotation) = paste("Test", 1:10, sep = "")

pheatmap(test, annotation = annotation)

 

Wednesday, January 30, 2013

R negative binomial distribution




> par(mfrow=c(5,2))
> for (i in (1:10)/10) {
+   barplot(dbinom(0:10, 10, i), main=sprintf('prob=%s',i))
+ }