Load packages
library(mosaic)
library(dplyr)
library(ggplot2)
library(ggvis)
library(parallel)


Simple Linear Regression

## Load cat dataset
library(MASS)
data(cats)
cats <- cats

## Scatterplot
cats %>% 
  ggvis(x=~Bwt, y=~Hwt, fill := "steelblue", stroke:="steelblue", fillOpacity:=.4, strokeOpacity:=.4) %>%
  layer_points() %>%
  add_axis("x", title = "Body weight (kg)") %>%
  scale_numeric("x", domain=c(1.5, 4), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Heart weight (g)") %>%
  scale_numeric("y", domain=c(6, 22), nice=FALSE, clamp=TRUE)

## Correlation
cor(cats$Bwt, cats$Hwt)
## [1] 0.8041274
## summary
cats %>%
  summarize(mean_bwt=mean(Bwt), sd_bwt=sd(Bwt), 
            mean_hwt=mean(Hwt), sd_hwt=sd(Hwt))
##   mean_bwt    sd_bwt mean_hwt   sd_hwt
## 1 2.723611 0.4853066 10.63056 2.434636
## Scatterplot with line
cats %>% 
  ggvis(x=~Bwt, y=~Hwt, fill := "steelblue", stroke:="steelblue", fillOpacity:=.4, strokeOpacity:=.4) %>%
  layer_points() %>%
  layer_model_predictions(model = "lm", strokeWidth:=5, stroke:="red") %>%
  add_axis("x", title = "Body weight (kg)") %>%
  scale_numeric("x", domain=c(1.5, 4), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Heart weight (g)") %>%
  scale_numeric("y", domain=c(6, 22), nice=FALSE, clamp=TRUE)

## Get coefficients for best-fitting line
lm(Hwt~Bwt, data=cats)
## 
## Call:
## lm(formula = Hwt ~ Bwt, data = cats)
## 
## Coefficients:
## (Intercept)          Bwt  
##     -0.3567       4.0341
## Load HSGPA SAU GPA data 
actgpa <- read.csv("http://www.bradthiessen.com/html5/data/f13s14.csv")

## Select ACT and GPA data
actgpa <-
  actgpa %>%
  dplyr::select(S1_P49HSGPA, S1_P49FallTermGPA, S1_P49ACTComposite, 
         S1_P49ACTMath, S1_P49ACTReading)
actgpa2 <- actgpa[complete.cases(actgpa),]

## Scatterplot
actgpa2 %>% 
  ggvis(x=~S1_P49HSGPA, y=~S1_P49FallTermGPA, fill := "steelblue", stroke:="steelblue", fillOpacity:=.2, strokeOpacity:=.4) %>%
  layer_points() %>%
  add_axis("x", title = "High School GPA") %>%
  scale_numeric("x", domain=c(1.8, 4), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "SAU 1st semester GPA") %>%
  scale_numeric("y", domain=c(0, 4.1), nice=FALSE, clamp=TRUE)

## Correlation and summary statistics
cor(S1_P49FallTermGPA, S1_P49HSGPA, data=actgpa2)
## [1] 0.6784345
mean(S1_P49FallTermGPA, data=actgpa2)
## [1] 2.63624
sd(S1_P49FallTermGPA, data=actgpa2)
## [1] 0.750968
mean(S1_P49HSGPA, data=actgpa2)
## [1] 3.235433
sd(S1_P49HSGPA, data=actgpa2)
## [1] 0.5436225
## Scatterplot with linear prediction
actgpa2 %>% 
  ggvis(x=~S1_P49HSGPA, y=~S1_P49FallTermGPA, fill := "steelblue", stroke:="steelblue", fillOpacity:=.2, strokeOpacity:=.4) %>%
  layer_points() %>%
    layer_model_predictions(model = "lm", strokeWidth:=5, stroke:="red") %>%
  add_axis("x", title = "High School GPA") %>%
  scale_numeric("x", domain=c(1.8, 4), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "SAU 1st semester GPA") %>%
  scale_numeric("y", domain=c(0, 4.1), nice=FALSE, clamp=TRUE)

## Coefficients of linear model
lm(S1_P49FallTermGPA~S1_P49HSGPA, data=actgpa2)
## 
## Call:
## lm(formula = S1_P49FallTermGPA ~ S1_P49HSGPA, data = actgpa2)
## 
## Coefficients:
## (Intercept)  S1_P49HSGPA  
##     -0.3960       0.9372
## Input beer dataset
brand <- c("Busch", "MGD", "Bud Light", "Coors Light", "Miller Lite", "Budweiser")
media <- c(8.7, 21.5, 68.7, 76.6, 100.1, 120.0)
ship <- c(8.1, 5.6, 20.7, 13.2, 15.9, 36.3)
beer <- data.frame(brand, media, ship)

## Scatterplot
beer %>% 
  ggvis(x=~media, y=~ship, fill := "steelblue", stroke:="steelblue", fillOpacity:=.9, strokeOpacity:=1) %>%
  layer_points(size:=100) %>%
  add_axis("x", title = "Media expenditures (millions of $)") %>%
  scale_numeric("x", domain=c(0, 125), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Shipments (millions of bottles)") %>%
  scale_numeric("y", domain=c(0, 40), nice=FALSE, clamp=TRUE)

## Scatterplot with least squares line
beer %>% 
  ggvis(x=~media, y=~ship, fill := "steelblue", stroke:="steelblue", fillOpacity:=.9, strokeOpacity:=1) %>%
  layer_points(size:=100) %>%
  layer_model_predictions(model = "lm", strokeWidth:=3, stroke:="red") %>%
  add_axis("x", title = "Media expenditures (millions of $)", grid=F) %>%
  scale_numeric("x", domain=c(0, 125), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Shipments (millions of bottles)", grid=F) %>%
  scale_numeric("y", domain=c(0, 40), nice=FALSE, clamp=TRUE)

## Scatterplot with lowess regression (curve)
beer %>% 
  ggvis(x=~media, y=~ship, fill := "steelblue", stroke:="steelblue", fillOpacity:=.9, strokeOpacity:=1) %>%
  layer_points(size:=100) %>%
  layer_smooths(strokeWidth:=5, stroke:="red") %>%
  add_axis("x", title = "Media expenditures (millions of $)") %>%
  scale_numeric("x", domain=c(0, 125), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Shipments (millions of bottles)") %>%
  scale_numeric("y", domain=c(0, 40), nice=FALSE, clamp=TRUE)

## Summary statistics for beer data
beer %>%
  summarize(avgShip = mean(ship), sdShip = sd(ship),
            acgMedia = mean(media), sdMedia = sd(media))
##    avgShip   sdShip acgMedia  sdMedia
## 1 16.63333 11.04711 65.93333 43.50166
cor(media, ship, data=beer)
## [1] 0.8288211
## Get coefficients for best-fitting line
lm(ship~media, data=beer)
## 
## Call:
## lm(formula = ship ~ media, data = beer)
## 
## Coefficients:
## (Intercept)        media  
##      2.7559       0.2105
## Create an example in which the best-fitting line
## doesn't fit well at all
x <- c(-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5)
y <- x^2
xy <- data.frame(x, y)

## Create scatterplot
xy %>% 
  ggvis(x=~x, y=~y, fill := "steelblue", stroke:="steelblue", fillOpacity:=.9, strokeOpacity:=1) %>%
  layer_points(size:=100) %>%
  layer_model_predictions(model = "lm", strokeWidth:=3, stroke:="red") %>%
  add_axis("x", title = "X", grid=F) %>%
  scale_numeric("x", domain=c(-5, 5), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Y", grid=F) %>%
  scale_numeric("y", domain=c(0, 25), nice=FALSE, clamp=TRUE)

## Create data with no correlation
x <- c(0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6)
y <- c(0, 1, 2, 3, 2, 1, 0, 0, -1, -2, -3, -2, -1, 0)
xy <- data.frame(x, y)

## Scatterplot
xy %>% 
  ggvis(x=~x, y=~y, fill := "steelblue", stroke:="steelblue", fillOpacity:=.9, strokeOpacity:=1) %>%
  layer_points(size:=100) %>%
  add_axis("x", title = "X", grid=F) %>%
  scale_numeric("x", domain=c(0, 6), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Y", grid=F) %>%
  scale_numeric("y", domain=c(-5, 5), nice=FALSE, clamp=TRUE)

# Best-fitting line
lm(y~x, data=xy)
## 
## Call:
## lm(formula = y ~ x, data = xy)
## 
## Coefficients:
## (Intercept)            x  
##           0            0
#########
## Get "full" output for beer regression
model = lm(ship~media, data=beer)
summary(model)
## 
## Call:
## lm(formula = ship ~ media, data = beer)
## 
## Residuals:
##      1      2      3      4      5      6 
##  3.513 -1.681  3.484 -5.678 -7.925  8.287 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2.75591    5.46812   0.504   0.6408  
## media        0.21048    0.07104   2.963   0.0414 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.911 on 4 degrees of freedom
## Multiple R-squared:  0.6869, Adjusted R-squared:  0.6087 
## F-statistic: 8.777 on 1 and 4 DF,  p-value: 0.04145
coef(model)
## (Intercept)       media 
##   2.7559140   0.2104765
confint(model)
##                    2.5 %     97.5 %
## (Intercept) -12.42603485 17.9378628
## media         0.01322851  0.4077246
vcov(model)
##             (Intercept)        media
## (Intercept)  29.9003908 -0.332776144
## media        -0.3327761  0.005047161
residuals(model)
##         1         2         3         4         5         6 
##  3.512940 -1.681159  3.484348 -5.678416 -7.924615  8.286902
fitted(model)
##         1         2         3         4         5         6 
##  4.587060  7.281159 17.215652 18.878416 23.824615 28.013098
deviance(model)
## [1] 191.0244
predict(model)
##         1         2         3         4         5         6 
##  4.587060  7.281159 17.215652 18.878416 23.824615 28.013098
anova(model)
## Analysis of Variance Table
## 
## Response: ship
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## media      1 419.17  419.17  8.7773 0.04145 *
## Residuals  4 191.02   47.76                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(model)
## [1] 43.79111
model.matrix(model)
##   (Intercept) media
## 1           1   8.7
## 2           1  21.5
## 3           1  68.7
## 4           1  76.6
## 5           1 100.1
## 6           1 120.0
## attr(,"assign")
## [1] 0 1
plot(model)

#########



### Randomization-based test of the slope of our regression line
randomizations <- do(10000) * coef(lm(S1_P49FallTermGPA~shuffle(S1_P49HSGPA), data=actgpa2))

histogram(~S1_P49HSGPA, data = randomizations, col="grey", xlim=c(-.4, 0.5), xlab = "Slope", width=.005)

tally(~ (S1_P49HSGPA >= 0.9372), format="prop", data=randomizations)
## 
##  TRUE FALSE 
##     0     1
####################
## Violin example ##
####################
violin <- read.csv("http://www.bradthiessen.com/html5/data/violin.csv")

## Scatterplot
violin %>% 
  ggvis(x=~years, y=~neural, fill := "steelblue", stroke:="steelblue", fillOpacity:=.8, strokeOpacity:=.8) %>%
  layer_points(size:=100) %>%
  add_axis("x", title = "Years played") %>%
  scale_numeric("x", domain=c(0, 20), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Neural activity") %>%
  scale_numeric("y", domain=c(0, 30), nice=FALSE, clamp=TRUE)

violin %>%
  summarize(mean=mean(years), sd=sd(years), mean2=mean(neural), sd2=sd(neural))
##   mean      sd    mean2      sd2
## 1  7.2 7.24273 15.56667 7.782459
cor(years, neural, data=violin)
## [1] 0.9279869
model <- lm(neural~years, data=violin)
summary(model)
## 
## Call:
## lm(formula = neural ~ years, data = violin)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8644 -2.3730  0.1614  2.3713  4.6471 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.3873     1.1149   7.523 4.35e-06 ***
## years         0.9971     0.1110   8.980 6.18e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.009 on 13 degrees of freedom
## Multiple R-squared:  0.8612, Adjusted R-squared:  0.8505 
## F-statistic: 80.63 on 1 and 13 DF,  p-value: 6.178e-07
anova(model)
## Analysis of Variance Table
## 
## Response: neural
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## years      1 730.21  730.21  80.633 6.178e-07 ***
## Residuals 13 117.73    9.06                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(model)

######################
## Prestige example ##
######################


## Load data
prestige <- read.csv("http://www.bradthiessen.com/html5/data/prestige.csv")

## Set occupation type as a factor variable
prestige$type <- as.factor(prestige$occ_type)

## Eliminate the job titles
prestige2 <- prestige %>%
  dplyr::select(-title, -blue, -professional, -white)

## Create scatterplot matrix with histograms on diagonals
panel.hist <- function(x, ...)
{
    usr <- par("usr"); on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts; y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
pairs(prestige2,
      cex = 1, pch = 1, bg = "steelblue",
      diag.panel = panel.hist, cex.labels = 1, font.labels = 2)

## ANOVA, summary statistics, Bartlett's test, Tukey/Bonferroni
model <- lm(prestige ~ as.factor(occ_type), data=prestige)
anova(model)
## Analysis of Variance Table
## 
## Response: prestige
##                     Df Sum Sq Mean Sq F value    Pr(>F)    
## as.factor(occ_type)  2  19467  9733.6  92.405 < 2.2e-16 ***
## Residuals           99  10428   105.3                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
prestige %>%
  group_by(occ_type) %>%
  summarize(n=n(), mean=mean(prestige), sd=sd(prestige))
## Source: local data frame [3 x 4]
## 
##   occ_type  n     mean        sd
## 1        0 49 36.08571 11.347320
## 2        1 23 42.24348  9.515816
## 3        2 30 67.90667  8.819255
bartlett.test(prestige ~ occ_type, data = prestige)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  prestige by occ_type
## Bartlett's K-squared = 2.4469, df = 2, p-value = 0.2942
TukeyHSD(model, conf.level=.95)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = x)
## 
## $`as.factor(occ_type)`
##          diff         lwr      upr     p adj
## 1-0  6.157764 -0.01491793 12.33045 0.0507001
## 2-0 31.820952 26.15954398 37.48236 0.0000000
## 2-1 25.663188 18.89483548 32.43154 0.0000000
plot(TukeyHSD(model))

with(prestige, pairwise.t.test(prestige, occ_type, p.adjust.method="bonferroni"))
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  prestige and occ_type 
## 
##   0       1      
## 1 0.059   -      
## 2 < 2e-16 4.4e-14
## 
## P value adjustment method: bonferroni
##### Construct linear model: Prestige = f(income)

## Scatterplot with regression line
prestige %>% 
  ggvis(x=~income, y=~prestige, fill := "steelblue", stroke:="steelblue", fillOpacity:=.3, strokeOpacity:=.5) %>%
  layer_points() %>%
    layer_model_predictions(model = "lm", strokeWidth:=5, stroke:="red") %>%
  add_axis("x", title = "Income") %>%
  scale_numeric("x", domain=c(0, 28000), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Prestige") %>%
  scale_numeric("y", domain=c(0, 100), nice=FALSE, clamp=TRUE)

## Create model
model <- lm(prestige~income, data=prestige)
# Get coefficients
coef(model)
##  (Intercept)       income 
## 27.141176368  0.002896799
# Summarize model to get R-squared
summary(model)
## 
## Call:
## lm(formula = prestige ~ income, data = prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.007  -8.378  -2.378   8.432  32.084 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.714e+01  2.268e+00   11.97   <2e-16 ***
## income      2.897e-03  2.833e-04   10.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.09 on 100 degrees of freedom
## Multiple R-squared:  0.5111, Adjusted R-squared:  0.5062 
## F-statistic: 104.5 on 1 and 100 DF,  p-value: < 2.2e-16
# Get ANOVA summary table
anova(model)
## Analysis of Variance Table
## 
## Response: prestige
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## income      1  15279 15279.3  104.54 < 2.2e-16 ***
## Residuals 100  14616   146.2                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Get confidence intervals for coefficients
confint(model)
##                    2.5 %       97.5 %
## (Intercept) 22.642116976 31.640235760
## income       0.002334692  0.003458907
### Create new data (job with income = 7000)
newdata = data.frame(income=7000)

## Find the predicted prestige for that new job
predict(model, newdata) 
##        1 
## 47.41877
## Find confidence interval
predict(model, newdata, interval="confidence") 
##        fit      lwr      upr
## 1 47.41877 45.04112 49.79642
## Find prediction interval
predict(model, newdata, interval="predict") 
##        fit      lwr      upr
## 1 47.41877 23.31552 71.52202
## Assumptions
plot(model)

## Homoscedasticity test
ncvTest(model)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 3.088455    Df = 1     p = 0.07884962
# Normality visualization: component + residual plot 
crPlots(model)

# Square transformation
lm(prestige^2 ~ income, data=prestige)
## 
## Call:
## lm(formula = prestige^2 ~ income, data = prestige)
## 
## Coefficients:
## (Intercept)       income  
##    447.7104       0.2999
## Robust regression
library(MASS)
rlm(prestige ~ income, data=prestige)
## Call:
## rlm(formula = prestige ~ income, data = prestige)
## Converged in 6 iterations
## 
## Coefficients:
##  (Intercept)       income 
## 25.121145855  0.003125411 
## 
## Degrees of freedom: 102 total; 100 residual
## Scale estimate: 13
## Quantile regression
library(quantreg)
rq(prestige ~ income, tau=.5, data=prestige)
## Call:
## rq(formula = prestige ~ income, tau = 0.5, data = prestige)
## 
## Coefficients:
##  (Intercept)       income 
## 23.945837640  0.003029339 
## 
## Degrees of freedom: 102 total; 100 residual
summary(rq(prestige ~ income, tau=.5, data=prestige))
## Warning in rq.fit.br(x, y, tau = tau, ci = TRUE, ...): Solution may be
## nonunique
## 
## Call: rq(formula = prestige ~ income, tau = 0.5, data = prestige)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 23.94584     15.10204 26.72000
## income       0.00303      0.00289  0.00453
## Lowess
lowess(prestige$prestige ~ prestige$income)
## $x
##   [1]   611   918  1656  1890  2370  2448  2594  2847  2901  3000  3016
##  [12]  3116  3148  3161  3472  3485  3485  3582  3617  3643  3739  3910
##  [23]  3930  3942  4036  4075  4199  4224  4330  4348  4443  4549  4614
##  [34]  4686  4696  4741  4753  4761  5052  5092  5134  5134  5180  5299
##  [45]  5449  5511  5562  5648  5795  5811  5902  5959  6112  6197  6259
##  [56]  6336  6462  6477  6565  6573  6686  6860  6928  6992  7059  7147
##  [67]  7405  7482  7562  7716  7869  7956  8034  8043  8049  8131  8206
##  [78]  8258  8316  8403  8425  8780  8845  8865  8880  8891  8895  9271
##  [89]  9593 10432 11023 11030 11377 12351 12480 14032 14163 14558 17498
## [100] 19263 25308 25879
## 
## $y
##   [1] 20.51871 21.82813 24.93925 25.91418 27.90203 28.22461 28.82842
##   [8] 29.87773 30.10282 30.51549 30.58219 31.00135 31.13548 31.18997
##  [15] 32.50728 32.56291 32.56291 32.97803 33.12781 33.23908 33.65130
##  [22] 34.38259 34.46812 34.51943 34.91100 35.07346 35.57983 35.68192
##  [29] 36.10845 36.18088 36.56315 36.97804 37.23245 37.51426 37.55561
##  [36] 37.74170 37.79132 37.82440 38.96850 39.12933 39.29821 39.29821
##  [43] 39.48317 39.96165 40.58139 40.83755 41.06556 41.45005 42.13674
##  [50] 42.21148 42.69183 42.99272 43.72096 44.12554 44.44584 44.84364
##  [57] 45.45230 45.52476 45.94987 45.98851 46.48047 47.21826 47.50658
##  [64] 47.75696 48.01906 48.36332 49.39793 49.72322 50.06118 50.72886
##  [71] 51.46755 51.88759 52.30803 52.35654 52.38888 52.83088 53.23514
##  [78] 53.51203 53.82086 54.28411 54.40125 56.40358 56.76581 56.87726
##  [85] 56.96085 57.02215 57.04444 59.08704 60.70490 63.38973 65.19523
##  [92] 65.21742 66.25486 68.37407 68.58814 71.42799 71.70646 72.56689
##  [99] 75.02568 76.11860 82.53114 83.12620
prestige %>% 
  ggvis(x=~income, y=~prestige, fill := "steelblue", stroke:="steelblue", fillOpacity:=.4, strokeOpacity:=.5) %>%
  layer_points() %>%
  layer_smooths(stroke:="red", strokeWidth:=5) %>%
  add_axis("x", title = "Income") %>%
  scale_numeric("x", domain=c(0, 28000), nice=FALSE, clamp=TRUE)  %>%
  add_axis("y", title = "Prestige") %>%
  scale_numeric("y", domain=c(0, 100), nice=FALSE, clamp=TRUE)

## Nonlinear -- quadratic model
quad = lm(prestige ~ income + I(income^2), data=prestige)
coef(quad)
##   (Intercept)        income   I(income^2) 
##  1.418319e+01  6.153506e-03 -1.433097e-07
summary(quad)
## 
## Call:
## lm(formula = prestige ~ income + I(income^2), data = prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.963  -7.967  -2.303   7.847  32.928 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.418e+01  3.515e+00   4.035 0.000108 ***
## income       6.154e-03  7.593e-04   8.104 1.43e-12 ***
## I(income^2) -1.433e-07  3.141e-08  -4.562 1.45e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.04 on 99 degrees of freedom
## Multiple R-squared:  0.596,  Adjusted R-squared:  0.5879 
## F-statistic: 73.03 on 2 and 99 DF,  p-value: < 2.2e-16
plot(prestige~income, data=prestige)
curve(14.18319+0.0061535*x-.0000001433097*x^2,add=T)

plot(quad)

## Log model
logmod = lm(prestige ~ log(income), data=prestige)
summary(logmod)
## 
## Call:
## lm(formula = prestige ~ log(income), data = prestige)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.114  -9.342  -1.218   8.101  30.454 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -139.856     16.954  -8.249  6.6e-13 ***
## log(income)   21.556      1.953  11.037  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.61 on 100 degrees of freedom
## Multiple R-squared:  0.5492, Adjusted R-squared:  0.5447 
## F-statistic: 121.8 on 1 and 100 DF,  p-value: < 2.2e-16
plot(prestige~income, data=prestige)
curve(-139.856 + 21.556*log(x),add=T)

plot(logmod)