Introduction

Libraries

Library that I used

library(ellipse)
library(ggplot2)
library(coda)
library(runjags)
library(rjags)
library(BayesFactor)
library(jagsUI)
library(HDInterval)  # for high density interval

Data

You can find the data on the course web-site. You need to change the path to where ever you put the data.

df <- read.table("C:/Users/cja/Dropbox/edps 590BAY/Lectures/7 Multiple Linear Regression/getting_what_paid_for_data.txt", header=TRUE)
names(df)
##  [1] "state"     "exp_pp"    "ave_pt"    "salary"    "taking"    "ave_v"    
##  [7] "ave_m"     "ave_tot"   "region"    "state_abv"
head(df)
##        state exp_pp ave_pt salary taking ave_v ave_m ave_tot region state_abv
## 1    Alabama  4.405   17.2 31.144      8   491   538    1029      S        AL
## 2     Alaska  8.963   17.6 47.951     47   445   489     934      W        AK
## 3    Arizona  4.778   19.3 32.175     27   448   496     944      W        AZ
## 4   Arkansas  4.459   17.1 28.934      6   482   523    1005      S        AR
## 5 California  4.992   24.0 41.078     45   417   485     902     CA        CA
## 6   Colorado  5.443   18.4 34.571     29   462   518     980      W        CO

Plot the Data

..and overlay the normal distribution. There are a number of ways to do this and what’s below keeps things in terms of frequencies. The first line sets up the histogram and saves it as “h”. The next two lines get data for horizontal and vertical axis. The third command puts density values of vertical fit values into scale of frequencies and uses information from the historgram. The “lines” command puts in the normal curve.

h =hist(df$ave_tot, breaks = 10, density = 10,
            col = "darkblue", 
            xlab = "Average Total SAT", 
            ylab = "Frequency",
            main="Normal curve Over histogram")
xfit <- seq(min(df$ave_tot), max(df$ave_tot), length = 50) 
yfit <- dnorm(xfit, mean = mean(df$ave_to), sd = sd(df$ave_tot)) 
yfit <- yfit * diff(h$mids[1:2]) * length(df$ave_tot) 
lines(xfit, yfit, col = "black", lwd = 2)

Sample Statistics

We need the sample size and observed sample mean. We will pretend that we know what the variance equals and use the sample variance for this.

n <- nrow(df)
ybar <- mean(df$ave_tot)
std <- sd(df$ave_tot)
s2 <- var(df$ave_tot)


descriptive <- c(n, ybar, std, s2)
names(descriptive) <- c("N","Ybar","std","s^2")
descriptive
##          N       Ybar        std        s^2 
##   50.00000  965.92000   74.82056 5598.11592
xs <- data.frame(df$ave_tot, df$exp_pp, df$salary, df$taking, df$ave_pt)
ctab <- cor(xs)
round(ctab,2)     # round correlations to 2 decimals for easier viewing
##            df.ave_tot df.exp_pp df.salary df.taking df.ave_pt
## df.ave_tot       1.00     -0.38     -0.44     -0.89      0.08
## df.exp_pp       -0.38      1.00      0.87      0.59     -0.37
## df.salary       -0.44      0.87      1.00      0.62      0.00
## df.taking       -0.89      0.59      0.62      1.00     -0.21
## df.ave_pt        0.08     -0.37      0.00     -0.21      1.00
# Figures of correlations in color 
colorfun <- colorRamp(c("#CC0000","white","#3366CC"), space="Lab")
plotcorr(ctab, col=rgb(colorfun((ctab+1)/2), maxColorValue=255),
         mar = c(0.1, 0.1, 0.1, 0.1))