library(tidyverse)

## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## filter(): dplyr, stats
## lag():    dplyr, stats

library(RCurl)

## Loading required package: bitops

## 
## Attaching package: 'RCurl'

## The following object is masked from 'package:tidyr':
## 
##     complete

lovett <- read_csv(getURL("https://raw.githubusercontent.com/chrischizinski/SNR_R_Group/master/data/ExperimentalDesignData/chpt2/lovett.csv"))

head(lovett)

## # A tibble: 6 × 5
##       STREAM ELEVATN   SO4 SO4MOD    CL
##        <chr>   <int> <dbl>  <dbl> <dbl>
## 1 Santa_Cruz     680  50.6   50.6  15.5
## 2    Colgate     628  55.4   55.4  16.4
## 3     Halsey     625  56.5   56.5  17.1
## 4 Batavia_Hi     663  57.5   57.5  16.8
## 5 Windham_Ri     616  58.3   58.3  18.3
## 6 Silver_Spr     451  63.0   63.0  15.7

lovett_sum<-lovett %>% 
      summarise(Mean_SO4 = mean(SO4),
                Mean_SO4mod = mean(SO4MOD),
                Mean_CL = mean(CL),
                Median_SO4 = median(SO4),
                Median_SO4mod = median(SO4MOD),
                Median_CL = median(CL),
                TrimMean_SO4 = mean(SO4,trim=0.05),
                TrimMean_SO4mod = mean(SO4MOD,trim=0.05),
                TrimMean_CL = mean(CL,trim=0.05),
                SD_SO4 = sd(SO4),
                SD_SO4mod = sd(SO4MOD),
                SD_CL = sd(CL),
                IQR_SO4 = IQR(SO4),
                IQR_SO4mod = IQR(SO4MOD),
                IQR_CL = IQR(CL),
                MAD_SO4 = mad(SO4),
                MAD_SO4mod = mad(SO4MOD),
                MAD_CL = mad(CL),
                SE_SO4 = sd(SO4)/sqrt(length(SO4)),
                SE_SO4mod = sd(SO4MOD)/sqrt(length(SO4MOD)),
                SE_CL = sd(CL)/sqrt(length(CL))) %>% 
      mutate(CILO_SO4 = Mean_SO4 - 1.96*SE_SO4,
             CIHI_SO4 = Mean_SO4 + 1.96*SE_SO4,
             CILO_SO4mod = Mean_SO4mod - 1.96*SE_SO4mod,
             CIHI_SO4mod = Mean_SO4mod + 1.96*SE_SO4mod,
             CILO_CL = Mean_CL - 1.96*SE_CL,
             CIHI_CL = Mean_CL + 1.96*SE_CL) %>% 
      gather(Mean, value, Mean_SO4:CIHI_CL) %>% 
      separate(Mean, c("Type", "Class"), sep ="_") %>% 
      spread(Class, value)

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

names(lovett)

## [1] "STREAM"  "ELEVATN" "SO4"     "SO4MOD"  "CL"

hub_vals<-apply(lovett[,c("SO4","SO4MOD","CL")], 2, function(x) huber(x)$mu)

names(hub_vals)<- c("SO4","SO4mod","CL")

hub_vals2<-data.frame(Type = "Huber", rbind(hub_vals))

lovett_sum <- rbind(lovett_sum, hub_vals2)

lovett_sum

## # A tibble: 10 × 4
##        Type        CL       SO4    SO4mod
## *     <chr>     <dbl>     <dbl>     <dbl>
## 1      CIHI 26.727466 63.568146 72.327934
## 2      CILO 18.954585 60.278008 58.077195
## 3       IQR  6.700000  8.000000  8.000000
## 4       MAD  5.782140  6.375180  6.375180
## 5      Mean 22.841026 61.923077 65.202564
## 6    Median 20.500000 62.100000 62.100000
## 7        SD 12.383068  5.241558 22.703020
## 8        SE  1.982878  0.839321  3.635393
## 9  TrimMean 21.594595 61.954054 61.954054
## 10    Huber 20.539791 61.954055 61.954055