.Rhistory

dat=read.table("http://www.biostat.jhsph.edu/bstcourse/bio751/data/babies.data",header=TRUE)
set.seed(0)
smokers=sample(dat$bwt[dat$smoke==1],10)
nonsmokers=sample(dat$bwt[dat$smoke==0],10)
library(rafalib)
mypar()
stripchart(list(smokers,nonsmokers),vertical=TRUE,method="jitter",col=c(1,2),pch=15,group.names=c("smokers","nonsmokers"))
cat("observed difference = ",mean(smokers)-mean(nonsmokers)," ounces")
(-6.1)*2
(-6.1)*-2
library(MASS)
ginv(A)
A<-matrix(c(1,2,4,2,1,1,3,1,2),nrow=3)
ginv(A)
library(MASS)
A<-matrix(c(1,2,4,2,1,1,3,1,2),nrow=3)
A
ginv(A)
1^-1
1/(1^-1)
A^-1
(A^-1)*A
A*(A^-1)
A^-1
ginv(ginv(A))
1/(A^-1)
A<-matrix(c(1,2,4,2,1,1,3,1,2),nrow=3)
A
ginv(A)
A^-1
C<-matix(c(4332),nrow=2)
C<-matrix(c(4332),nrow=2)
C
C<-matrix(c(4,3,3,2),nrow=2)
C
C^-1
(C)^-1
2^2
2^3
2^-3
2^-2
1^-2
2^-2
2^(-2)
2*-2
1/(2^2)
2^(-2)
ginv(A)
mns <- vector(length=1e5)
for(i in seq(1e5)) mns[i] <- mean(rexp(50))
mns2 <- vector(length=1e5)
for(i in seq(1e5)) mns2[i] <- mean(rexp(10))
plot(density(mns))
lines(density(mns2), lty=2)
rm(i)
data(maPooling)
library(dagdata)
e <- data(maPooling)
rm(list=ls())
library(dagdata)
e <- data(maPooling)
rm(list=ls())
library(dagdata)
data(maPooling)
e <- maPooling
rm(list=ls())
e <- data(maPooling)
class(e)
dim(e)
dim(maPooling)
install.packages("data.table")
library(data.table)
rm(list=ls())
data(maPooling)
rm(list=ls())
e <- data.table(data(maPooling))
rm(list=ls())
data(maPooling)
e -data.table(maPooling)
e <- data.table(maPooling)
dim(e)
rm(list=ls())
e <- data.table(data(maPooling))
e <- data.table(as.data.frame(data(maPooling)))
rm(list=ls())
e <- data.table(as.data.frame(data(maPooling)))
e <- data.table(as.data.frame(maPooling))
rm(list=ls())
e <- data.table(as.data.frame(maPooling))
library(dagdata)
e <- data.table(as.data.frame(maPooling))
data.table(maPooling)
data(maPooling)
e <- data.table(as.data.frame(maPooling))
rm(list=ls())
data(maPooling)
e <- maPooling
dim(e)
e <- data.table(maPooling)
e <- data.table(as.data.frame(maPooling))
class(e)
e <- as.data.table(maPooling)
rm(list=ls())
data(maPooling)
e <- as.data.table(maPooling)
e <- data.table(as.data.frame(maPooling))
rm(list=ls())
e <- maPooling
data(maPooling)
e <- maPooling
e[1,1]
head(e)
library(data.table)
> data(maPooling)
rm(list=ls())
data(maPooling)
e <- data.table(data(maPooling))
e <- as.data.frame(data(maPooling))
e <- as.data.frame(data(maPooling))
e <- as.data.frame(maPooling)
head(e)
# Background Assessment
rm(list=ls())
library(dagdata)
data(maPooling)
u = exprs(maPooling)[,1]
v = exprs(maPooling)[,2]
x = exprs(maPooling)[,3]
y = exprs(maPooling)[,4]
head(u)
# names(u)
dim(u)
length(u)
class(u)
u = exprs(maPooling)[,1]
v = exprs(maPooling)[,2]
x = exprs(maPooling)[,3]
y = exprs(maPooling)[,4]
head(u)
# names(u)
dim(u)
length(u)
class(u)
validIP <- "((25[0-5]|2[0-4][0-9]|[01]?[0-9][0-9]?)[.]){3}(25[0-5]|2[0-4][0-9]|[01]?[0-9][0-9]?)"
any(grep(validIP, ipmessage)
)
system("ipconfig", intern = TRUE)
system("ipconfig", intern = TRUE)
rm(validIP)
a <- system("ipconfig", intern = TRUE)
library(sos)
# Let's find in the help pages for a function that includes the keywords
# folder AND path
a <- findFn("split")
# Let's take a look at the packages returned
summary(a, minPackages = 87)
# Now let's subset the results to extract only those functions contained at
# the package "base" which is the one we are working with
b <- a[a$Package=="base",]
# Open the filtered HTML page
b
d <- strsplit(a,":")
c <- system("ipconfig", intern = TRUE)
d <- strsplit(c,":")
d
a <- system("ipconfig", intern = TRUE)
a
d <- strsplit(c,": ")
d
d <- unlist(strsplit(c,": "))
d
c <- system("ipconfig", intern = TRUE)
d <- unlist(strsplit(c,": "))
d
c
d <- unlist(strsplit(c,": ")[2])
d
a <- findFn("internet connection")
# Let's take a look at the packages returned
summary(a, minPackages = 87)
# Now let's subset the results to extract only those functions contained at
# the package "base" which is the one we are working with
b <- a[a$Package=="base",]
# Open the filtered HTML page
b
a <- findFn("internet")
# Let's take a look at the packages returned
summary(a, minPackages = 87)
# Now let's subset the results to extract only those functions contained at
# the package "base" which is the one we are working with
b <- a[a$Package=="base",]
# Open the filtered HTML page
b
internet.info
internet.info()
a <- findFn("ipconfig")
a <- findFn("ip")
# Let's take a look at the packages returned
summary(a, minPackages = 87)
# Now let's subset the results to extract only those functions contained at
# the package "base" which is the one we are working with
b <- a[a$Package=="base",]
# Open the filtered HTML page
b
a <- findFn("test connection")
# Let's take a look at the packages returned
summary(a, minPackages = 87)
b <- a[a$Package=="base",]
# Open the filtered HTML page
b
showConnections()
showConnections(all=T)
url("www.google.com")
Sys.info()
system("ping www.google.com")
a <- findFn("ping")
# Let's take a look at the packages returned
summary(a, minPackages = 87)
# Now let's subset the results to extract only those functions contained at
# the package "base" which is the one we are working with
b <- a[a$Package=="base",]
# Open the filtered HTML page
b
?help.search
library(Biobase)
library(dagdata)
data(maPooling)
class(maPooling)
e <- maPooling
head(pData(e))
individuals <- which(rowSums(pData(e)) == 1)
head(maPooling)
library(dagdata);
data(maPooling);
u = exprs(maPooling)[,1];
v = exprs(maPooling)[,2];
x = exprs(maPooling)[,3];
y = exprs(maPooling)[,4] )
y = exprs(maPooling)[,4];
help.search("atomic vector", package="base")
help.search("atomic vectors", package="base")
help.search("atomic", package="base")
help.search("vector", package="base")
help.search("atomic", package="base")
help.search("lexical")
help.search("lexical", package="base")
base::scoping
?base::scoping
example(base::scoping)
example("base::scoping")
example(scoping)
example("scoping")
help.search("code demonstrations", package="base")
help.search("demonstrations", package="base")
?rownames
help.search("lexical",package="base")
demo(scoping)
cran.help <- findFn("{hadoop}",sortby="MaxScore");
library(sos)
cran.help <- findFn("{hadoop}",sortby="MaxScore");
View(cran.help)
disease  <- c(rep(0, 180), rep(1, 20), rep(0, 40), rep(1, 10))
genotype <- c(rep("AA", 200), rep("aa", 50))
tab      <- table(genotype, disease)
tab
is(tab)
names(tab)
rownames(tab)
browser(table)
browser(tab)
browser(table)
library(Rdocumentation)
?table
tab      <- as.data.frame(table(genotype, disease))
tab
tab      <- table(genotype, disease)
order(tab,c("AA","aa"))
tab
disease  <- c(rep(1, 20), rep(0, 180), rep(0, 40), rep(1, 10))
genotype <- c(rep("AA", 200), rep("aa", 50))
tab      <- table(genotype, disease)
tab
Biases, systematic errors and unexpected variability are common in genomics data. Failure to discover these problems often leads to flawed analyses and false discoveries. As an example, consider that experiments sometimes fail and not all data processing pipelines are designed to detect these. Yet, these pipelines still give you an answer and the from the final results it may be hard or impossible to notice an error was made. In later modules we will cover many other examples.
opts_chunk$set(fig.path="figure/exploratory_data_analysis-")
=======
opts_chunk$set(fig.path=paste0("figure/", sub("(.*).Rmd","\\1",basename(knitr:::knit_concord$get('infile'))), "-"))
library(UsingR)
x=father.son$fheight
opts_chunk$set(fig.path="figure/exploratory_data_analysis-")
Graphing data is a powerful approach to detecting these problems. We refer to this as _exploratory data analyis_ (EDA). Many important methodological contributions to genomics data analysis were initiated as discovery made via EDA. We will show some useful exploratory plots for gene expression data measured with microarrays and NGS. We start with a general introduction to EDA using height data.
Biases, systematic errors and unexpected variability are common in genomics data. Failure to discover these problems often leads to flawed analyses and false discoveries. As an example, consider that experiments sometimes fail and not all data processing pipelines are designed to detect these. Yet, these pipelines still give you an answer and the from the final results it may be hard or impossible to notice an error was made. In later modules we will cover many other examples.
library(UsingR)
x = father.son$fheight
round(sample(x, 20), 1)
round(sample(x, 20), 1)
round(sample(x, 20), 1)
round(sample(x, 20), 1)
round(sample(x, 20), 1)
round(sample(x, 20), 1)
hist(x, breaks = seq(floor(min(x)), ceiling(max(x))), main = "", xlab = "Height")
xs <- seq(floor(min(x)), ceiling(max(x)), 0.1)
plot(xs, ecdf(x)(xs), type = "l", xlab = "x=Height", ylab = "F(x)")
1 - pnorm(72, mean(x), sd(x))
ps <- seq(0.01, 0.99, 0.01)
qs <- quantile(x, ps)
normalqs <- qnorm(ps, mean(x), sd(x))
plot(normalqs, qs, xlab = "Normal percentiles", ylab = "Height percentiles")
abline(0, 1)  ##identity line
qqnorm(x)
qqline(x)
hist(exec.pay)
hist(exec.pay)
qqnorm(exec.pay)
qqline(exec.pay)
boxplot(exec.pay, ylab = "10,000s of dollars", ylim = c(0, 400))
hist(exec.pay)
qqnorm(exec.pay)
qqline(exec.pay)
boxplot(exec.pay, ylab = "10,000s of dollars", ylim = c(0, 400))
rm(list=ls())
data("father.son")
x = father.son$fheight
y = father.son$sheight
plot(x,
y,
xlab = "Father's height in inches",
ylab = "Son's height in inches",
main = paste("correlation =", signif(cor(x, y), 2)))
boxplot(split(y, round(x)))
rm(list=ls())
data("father.son")
x = father.son$fheight
y = father.son$sheight
plot(x,
y,
xlab = "Father's height in inches",
ylab = "Son's height in inches",
main = paste("correlation =", signif(cor(x, y), 2)))
boxplot(split(y, round(x)))
print(mean(y[round(x) == 72]))
x = (x - mean(x))/sd(x)
y = (y - mean(y))/sd(y)
means = tapply(y, round(x * 4)/4, mean)
fatherheights = as.numeric(names(means))
plot(fatherheights,
means,
ylab = "average of strata of son heights",
ylim = range(fatherheights))
abline(0, cor(x, y))
data("father.son")
x = father.son$fheight
y = father.son$sheight
boxplot(split(y, round(x)))
print(mean(y[round(x) == 72]))
x = (x - mean(x))/sd(x)
y = (y - mean(y))/sd(y)
means = tapply(y, round(x * 4)/4, mean)
fatherheights = as.numeric(names(means))
plot(fatherheights,
means,
ylab = "average of strata of son heights",
ylim = range(fatherheights))
abline(0, cor(x, y))
rm(list=ls())
a = rnorm(100)
a[1] = 10
b = rnorm(100)
b[1] = 11
plot(a, b, main = paste("correlation =", signif(cor(a, b), 2)))
a = rnorm(100)
a[1] = 10
b = rnorm(100)
b[1] = 11
plot(a, b, main = paste("correlation =", signif(cor(a, b), 2)))
rm(list=ls())
a = rnorm(100)
a[1] = 10
b = rnorm(100)
b[1] = 11
plot(a, b, main = paste("correlation =", signif(cor(a, b), 2)))
mylocal.path <- "C:/Users/Diego/Documents/GitHub/R-Notes";
setwd(mylocal.path)
dir
dir()
purl("Minitutorial_01.Rmd")
library(MASS)
purl("Minitutorial_01.Rmd")
dir()
library(MASS)
purl("Minitutorial_01.Rmd")
??purl
library(knitr)
purl("Minitutorial_01.Rmd")
dir()
source("Minitutorial_01.R")
source("Minitutorial_01.R")
help(Control)
library(help=base)
help(Control)
help(Control)
?Control
help("Control")
detach("package:Rdocumentation", unload=TRUE)
help("Control")
purl("Minitutorial_01.Rmd")
dir()
source("Minitutorial_01.R")
load("Minitutorial_01.R")
vignette("Sweave");
source("Minitutorial_01.R")
library(help=base)
help("Control")
apropos("help")
?help.search;
## ----localhelpex---------------------------------------------------------
help.search("remove objects", package = "base");
## ----ls,tidy=FALSE-------------------------------------------------------
?ls();
## ----find----------------------------------------------------------------
find("rm")
## ----examples------------------------------------------------------------
example(rm);
## ----browsevignettes-----------------------------------------------------
browseVignettes(all=T);
## ----vignette------------------------------------------------------------
vignette("Sweave");
## ----vigbrowser----------------------------------------------------------
browseVignettes("knitr");
## ----rdocumentation------------------------------------------------------
# Uncomment if you want to go to the web page
browseURL("http://www.rdocumentation.org/")
## ----rdocumentationgit---------------------------------------------------
# Uncomment if you want to go to the web page
browseURL("https://github.com/Data-Camp/Rdocumentation")
## ----devtools------------------------------------------------------------
# library(devtools)
# install_github("Rdocumentation","Data-Camp")
## ----RdocEx, tidy=FALSE--------------------------------------------------
library(Rdocumentation)
?rm
## ------------------------------------------------------------------------
detach("package:Rdocumentation", unload=TRUE)
purl("Minitutorial_02.Rmd")
dir()
purl("Minitutorial_03.Rmd")
dir()
methods(class="lm")
library(geneplotter)
geneplotter::
methods(class="data.frame")
?by.data.frame
?split.data.frame
g <- airquality$Month
l <- split(airquality, g)
head(g)
g
l
airquality
l
l <- lapply(l, transform, Oz.Z = scale(Ozone))
aq2 <- unsplit(l, g)
head(aq2)
with(aq2, tapply(Oz.Z,  Month, sd, na.rm = TRUE))
methods(class="data.frame")
?Ops
?Ops.data.frame
?transform.data.frame
library(sos); # Review the use of the require() function.
cran.help <- findFn("{file extensions}",sortby="MaxScore");
is(cran.help)
class(cran.help)
?is
is(cran.help, "data.frame")
?Summary
?Summary.data.frame
methods("findFn")
h <-summary(cran.help)
class(h)
?summary.findFn
h <-summary(cran.help, "findFn")
methods(cran.help, "findFn")
methods("findFn","sos")
?methods
showMethods("findFn")
getMethod("findFn","data.frame")
?summary.findFn
methods('summary.findFn', 'list')
?S3
S3(cran.help)
S3Part(cran.help)
S3Part(findFn)
S3Class(findFn)
getS3method(findFn)
getAnywhere(findFn)
showMethods(findFn)
getS3method("findFn","sos")
methods("findFn")
Methods("findFn")
testInheritedMethods("findFn")
testInheritedMethods(findFn)
testInheritedMethods(cran.help)
getGeneric(findFn)
getGeneric("findFn")
getMethodsMetaData("findFn")
methods(cran.help, findFn)
methods(findFn, findFn)
methods("findFn", "findFn")
help(showMethods)
showMethods(findFn)
showMethods(findFn, "data.frame")
sos::?
?.sos