## @knitr unnamed-chunk-1
(5^2)*(10-8)/3 + 1       ## An expression
x <- (5^2)*(10-8)/3 + 1  ## An assignment
x


## @knitr unnamed-chunk-2
x+1
n <- 50
x/n


## @knitr unnamed-chunk-3
class(x)


## @knitr unnamed-chunk-4
x <- 1:9  ## Creates a vector of numbers 1, 2, ..., 9
mean(x)
median(x)
sd(x)
min(x)
sum(x)
x^2
sum(x^2)


## @knitr unnamed-chunk-5
x <- rnorm(1000)
y <- rnorm(1000, mean=3)
hist(c(x, y))
hist(c(x, y), col="gray", border="white", breaks=40)


## @knitr unnamed-chunk-6
x1 <- c(-0.1, 0.3, -0.9, 0.2, -0.6, -2.3)
x2 <- c("red", "blue", "green", "purple")
x3 <- x1 > 0
x3


## @knitr unnamed-chunk-7
names(x1) <- c("Justine", "Rachel", "Meng", "Xiuhua", "Whitney", "Paul")
x1
names(x1)
x1["Justine"]


## @knitr unnamed-chunk-8
X <- matrix(x1, nrow=2)
X


## @knitr unnamed-chunk-9
matrix(x1, nrow=2, byrow=TRUE)


## @knitr unnamed-chunk-10
X1 <- matrix(rnorm(6), ncol=3)
X2 <- matrix(rnorm(6), ncol=3)
rbind(X1, X2)
cbind(X1, X2)


## @knitr unnamed-chunk-11
tt <- t.test(1:7, 5:10)
tt
tt$conf.int
tt$p.value


## @knitr unnamed-chunk-12
tips <- read.delim("http://web.as.uky.edu/statistics/users/pbreheny/760/data/tips.txt")
head(tips)
class(tips$TotBill)
class(tips$Sex)


## @knitr unnamed-chunk-13
with(tips, mean(TotBill))
mean(TotBill)
attach(tips)
mean(TotBill)
detach(tips)


## @knitr unnamed-chunk-14
tips$Rate <- with(tips, Tip/TotBill)


## @knitr unnamed-chunk-15
table(tips$Sex)
levels(tips$Sex)
barplot(table(tips$Day))


## @knitr unnamed-chunk-16
x1[x1 > 0]  ## Logical vector
x1[1:3]  ## Positive integers
x1[-(5:6)]  ## Negative integers
x1[c("Whitney", "Meng")]  ## Names
x1[]


## @knitr unnamed-chunk-17
tips[tips$Tip >= 7, ]


## @knitr unnamed-chunk-18
x <- sample(1:10, 3)
y <- sample(1:10, 6)
x
x + 2  ## Adds 2 to each element of x
x + x  ## Adds x to itself elementwise
x + y  ## Adds x to y elementwise, 'recycling' x


## @knitr unnamed-chunk-19
x <- TRUE
y <- FALSE
x & y
x | y


## @knitr unnamed-chunk-20
seq(0, 1, .2)
seq(0, 1, len=5)
rep(1, 5)
rep(1:2, 2)
rep(1:2, c(2,2))


## @knitr unnamed-chunk-21
sum(is.na(airquality$Ozone))
mean(airquality$Ozone)
mean(airquality$Ozone, na.rm=TRUE)


## @knitr unnamed-chunk-22
1/0
0/0


## @knitr unnamed-chunk-23
A <- matrix(rpois(12, 3), nrow=4)
B <- matrix(rpois(12, 3), nrow=4)
A*B + 2*A


## @knitr unnamed-chunk-24
t(A) %*% B
crossprod(A, B)
A %*% B


## @knitr unnamed-chunk-25
X <- crossprod(A)
X.inv <- solve(X)
X %*% X.inv


## @knitr unnamed-chunk-26
paste(LETTERS[1:4], 1:4)
paste(LETTERS[1:4], 1:4, sep="")
a <- paste(LETTERS[1:4], 1:4, sep="-")
nchar(a)
substr(a, 1, 2)


## @knitr unnamed-chunk-27
apply(airquality, 2, mean)
apply(airquality, 2, mean, na.rm=TRUE)


## @knitr unnamed-chunk-28
sample(1:10, 5)
sample(LETTERS[1:10], 5)
sample(LETTERS[1:5], 10, replace=TRUE)


## @knitr unnamed-chunk-29
ols <- function(XX, yy) {
  missing.data <- apply(is.na(XX), 1, any) | is.na(yy)
  X <- cbind(Intercept=1, as.matrix(XX[!missing.data,]))
  y <- yy[!missing.data]
  solve(crossprod(X)) %*% t(X) %*% y
}
ols(airquality[,-1], airquality$Ozone)
## Check against existing R function:
lm(Ozone~., airquality)


## @knitr unnamed-chunk-30
require(lattice)
xyplot(Tip~TotBill|Smoker, data=tips)


## @knitr unnamed-chunk-31
x <- airquality$Ozone
x <- tips$Sex
if (is.numeric(x)) sqrt(x) else stop("Can't take sqrt of things that aren't numbers")


## @knitr unnamed-chunk-32
total <- 0
for (i in 1:10) {
  total <- total + i
}
total


## @knitr unnamed-chunk-33
N <- 1000
n <- 10
SD <- 1:10
pW <- pS <- matrix(NA, nrow=N, ncol=length(SD))
for (i in 1:length(SD)) {
  print(i)
  for (j in 1:N) {
    s1 <- rnorm(n)
    s2 <- rnorm(n, sd=SD[i])
    pW[j,i] <- t.test(s1, s2, var.equal=FALSE)$p.value
    pS[j,i] <- t.test(s1, s2, var.equal=TRUE)$p.value
  }
}
plot(SD,apply(pW < .05, 2, mean), type="l", lwd=3, col="red", ylab="Type I error rate",xlab="Ratio of standard deviations", ylim=c(0,0.1))
lines(SD, apply(pS < .05, 2, mean), lwd=3, col="blue")
legend("topleft", col=c("red", "blue"), lwd=3, legend=c("Welch","Student"))
abline(h=.05, col="gray")