Skip to content

Instantly share code, notes, and snippets.

@marutter
Created October 7, 2010 18:59
Show Gist options
  • Save marutter/615679 to your computer and use it in GitHub Desktop.
Save marutter/615679 to your computer and use it in GitHub Desktop.
data <- read.csv("sparrow.csv")
attach(data)
plot(ppha,pairs,xlab="Pedestrians per ha per min",ylab="Breeding pairs per ha")
res <- lm(pairs~ppha)
summary(res)
abline(res)
# Add the quadratic term
res2 <- lm(pairs~ppha+ppha^2)
summary(res2)
res2 <- lm(pairs~ppha+I(ppha^2))
summary(res2)
abline(res2)
# Plot Results
nd <- data.frame(ppha=seq(min(ppha),max(ppha),(max(ppha)-min(ppha))/45))
plot(ppha,pairs,xlab="Pedestrians per ha per min",ylab="Breeding pairs per ha")
# Cubic
res3 <- lm(pairs~poly(ppha,3))
summary(res3)
matlines(nd,predict(res3,nd),type="l",col="red",lwd=2)
# Quartic
res4 <- lm(pairs~poly(ppha,4))
summary(res4)
matlines(nd,predict(res4,nd),type="l",col="green",lwd=2)
#
# Adjusted R-squared
#
summary(res3)$r.squared
summary(res4)$r.squared
summary(res3)$adj.r.squared
summary(res4)$adj.r.squared
#
# Higher Degree
#
nd <- data.frame(ppha=seq(min(ppha),max(ppha),(max(ppha)-min(ppha))/250))
resh <- lm(pairs~poly(ppha,10))
plot(ppha,pairs,xlab="Pedestrians per ha per min",ylab="Breeding pairs per ha")
matlines(nd,predict(resh,nd),type="l",col="red",lwd=2)
resh <- lm(pairs~poly(ppha,22))
plot(ppha,pairs,xlab="Pedestrians per ha per min",ylab="Breeding pairs per ha")
matlines(nd,predict(resh,nd),type="l",col="red",lwd=2)
summary(resh)$r.squared
summary(resh)$adj.r.squared
#
# What Polynomial?
#
res2 <- lm(pairs~ppha+I(ppha^2))
summary(res2)
res2b <- lm(pairs~poly(ppha,2))
summary(res2b)
res2c <- lm(pairs~poly(ppha,2,raw=TRUE))
summary(res2c)
res10 <- lm(pairs~poly(ppha,10,raw=TRUE))
summary(res10)
resorth <- lm(pairs~poly(ppha,10))
summary(resorth)
model.matrix(res2)
model.matrix(res2b)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment