Skip to content

Instantly share code, notes, and snippets.

@dwinter
Created February 28, 2011 11:59
Show Gist options
  • Save dwinter/847229 to your computer and use it in GitHub Desktop.
Save dwinter/847229 to your computer and use it in GitHub Desktop.
plotting and testing some of Ken Ring's theories
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment