plotting and testing some of Ken Ring's theories

gistfile1.txt

ring.r

library(ggplot)
library(arm)

eq <- read.csv('eq2.csv')
attach(eq)
eq$date <- as.Date(eq$date, "%d/%m/%Y")

#let's get plotting...

theme_set(theme_bw())
molten.eq <- melt(eq, id.var='date')
head(molten.eq)

#         date variable value
# 1 2010-09-04   lunar     4
# 2 2010-09-05   lunar     3
# 3 2010-09-06   lunar     2
# 4 2010-09-07   lunar     1
# 5 2010-09-08   lunar     0
# 6 2010-09-09   lunar     1

p <- ggplot(molten.eq, aes(date, value))
p + geom_line() + facet_grid(variable ~ ., scales="free")

#Not bad, but I want to smooth the moon series, I'm sure that's possible in ggplot
#but I just made two graphs and pasted them together in inkscape:

last_plot() + stat_smooth()

#now, how to Ring's predictors go?

pMoon <- ggplot(eq, aes(lunar, energy))
pMoon + geom_point() + stat_smooth(method="lm", formula= y ~ poly(x,2))

pOrbit <- ggplot(eq, aes(ld, energy))
pOrbit + geom_point() + stat_smooth(method="lm")
  
# lets test them then...
modKR <- glm(energy ~ I(date -as.Date("2010-09-04")) + ld + poly(lunar, 2) + ld:poly(lunar,2))

# Then a bunch of model selection with AIC to get
modKR <- glm(energy ~ ld + I(date -as.Date("2010-09-04")))


#how about a predictive model?
#overdispresion parameter on full model with quasibinomial was 1.1, so didn't bother with it...

modPredict <- glm(large.10 ~ I(date -as.Date("2010-09-04")) + ld + poly(lunar, 2) + ld:poly(lunar,2), 
                  family="binomial", na.action=na.exclude)

#... and after model selection
modPredict <- glm(large.10 ~ I(date -as.Date("2010-09-04")), family="binomial", na.action=na.exclude)