# STAT3530 -- getting started in R
# Run this a line at a time and look at what each command gives you.

# --- load a dataset -------------------------------------------------------
# this reads the data from the course website and creates an object
# called `prevend`
load(url("https://tdruiz-stat3530.share.connect.posit.cloud/_data/prevend.RData"))

head(prevend)          # first few rows
nrow(prevend)          # number of observations

# --- numerical summaries --------------------------------------------------
mean(prevend$rfft)     # sample mean
var(prevend$rfft)      # sample variance
sd(prevend$rfft)       # sample standard deviation

# --- one variable at a time -----------------------------------------------
hist(prevend$rfft)

# --- two variables at a time ----------------------------------------------
plot(prevend$age, prevend$rfft)

# --- fit a simple linear regression ---------------------------------------
fit <- lm(rfft ~ age, data = prevend)   # response ~ predictor

fit                    # the least squares estimates
sigma(fit)             # the estimate of sigma

abline(fit, col = 'blue')     # add the fitted line to the scatterplot

# --- your turn ------------------------------------------------------------
# `cars` is built into R: the speed of 50 cars and the distance each took
# to stop
head(cars)

# using only the commands above, and the cars data:
#   1. make a scatterplot of dist against speed
#   2. fit the model dist ~ speed with lm()
#   3. find the estimates of beta0, beta1, and sigma
#   4. add the fitted line to your scatterplot

# --- vectors and matrices -------------------------------------------------
# the same a, b, A, B we used in the slides
a <- c(-1, 0, 2)
b <- c(3, 5, 6)

A <- matrix(c(1, 0, 3,
              4, -2, 0), nrow = 2, byrow = TRUE)
B <- matrix(c(0, 1,
              -1, 2,
              1, 0), nrow = 3, byrow = TRUE)

dim(A)                 # rows, columns
t(A)                   # transpose

# careful: * multiplies entry by entry, %*% is matrix multiplication
a * b                  # three products, no sum
sum(a * b)             # the inner product
t(a) %*% b             # the same number, as a 1x1 matrix

A %*% B                # 2x3 times 3x2 gives 2x2
B %*% A                # 3x2 times 2x3 gives 3x3
# A %*% A              # error: dimensions don't match

diag(3)                # the 3x3 identity
solve(A %*% B)         # solve() inverts a square matrix

# --- the slope estimate, by hand ------------------------------------------
# in the slides we wrote the slope estimate as a weighted sum of responses
x <- prevend$age
y <- prevend$rfft

w <- (x - mean(x)) / sum((x - mean(x))^2)

sum(w)                 # the weights sum to zero
beta1 <- sum(w * y)
beta0 <- mean(y) - beta1 * mean(x)

c(beta0, beta1)        # compare with
coef(fit)              # what lm() reported

# --- simulating from a normal distribution --------------------------------
set.seed(82226)        # so the random draws repeat next time

z <- rnorm(1000)       # 1000 draws from N(0, 1)

mean(z)                # close to 0, but not exactly
sd(z)                  # close to 1

hist(z)

# rnorm() also takes a mean and an sd; try changing them
scores <- rnorm(1000, mean = 100, sd = 15)
hist(scores)

# --- simulating from the simple linear model ------------------------------
# in matrix form the model is  Y = X * beta + e
n <- 50
x_sim <- 1:n           # predictors: we choose these, they are not random

X <- cbind(1, x_sim)   # design matrix: a column of ones, then x
beta <- c(10, 2)       # the true intercept and slope

head(X)                # first few rows
dim(X)                 # n rows, 2 columns

e <- rnorm(n, mean = 0, sd = 5)   # errors: this is the random part
y_sim <- X %*% beta + e           # matrix times vector, plus the errors

plot(x_sim, y_sim)
abline(10, 2, col = 'red')                # the truth
abline(lm(y_sim ~ x_sim), col = 'blue')   # our estimate

coef(lm(y_sim ~ x_sim))   # close to 10 and 2, but not equal

# run these three lines a few times: the estimate changes every time,
# because the response is random
e <- rnorm(n, mean = 0, sd = 5)
y_sim <- X %*% beta + e
coef(lm(y_sim ~ x_sim))