kd <- read.csv("http://www.stat.ufl.edu/~winner/data/kentuckyderby.csv", header=TRUE) attach(kd); names(kd) Year125 <- Year[Length==1.25] Time125 <- Time[Length==1.25] Length125 <- Length[Length==1.25] Year125.0 <- Year125-min(Year125) Speed125 <- 1609.34*Length125/Time125 summary(Speed125) Speed125.2008 <- Speed125[Year125<=2008] Year125.2008 <- Year125[Year125<=2008] kd.mod1 <- nls(Speed125 ~ b0 + (b1-b0)*exp(b2*(Year125.0-b3))/ (1+exp(b2*(Year125.0-b3))), start=c(b0=1,b1=20,b2=1,b3=45)) summary(kd.mod1) AIC(kd.mod1) confint(kd.mod1) plot(Year125,Speed125) lines(Year125,predict(kd.mod1,Year125)) plot(Year125,resid(kd.mod1)) qqnorm(resid(kd.mod1)); qqline(resid(kd.mod1)) kd.mod2 <- nls(Time125 ~ b0 + (b1*Year125.0)/(b2+Year125.0), start=c(b0=130,b1=-15,b2=45)) summary(kd.mod2) AIC(kd.mod2) plot(Year125.0,Time125) lines(Year125.0,predict(kd.mod2,Year125.0)) ### Matrix form for logistic model Y <- Speed125 X.t <- Year125.0 beta.old <- c(10,20,0.1,60) diff.beta <- 1000 while (diff.beta > .00001) { exp.t <- exp(beta.old[3] * (X.t - beta.old[4])) G1 <- 1 - (exp.t/(1+exp.t)) G2 <- exp.t/(1+exp.t) G3 <- (beta.old[2]-beta.old[1]) * (X.t-beta.old[4]) * (exp.t/(1+exp.t)^2) G4 <- (beta.old[2]-beta.old[1]) * (-beta.old[3]) * (exp.t/(1+exp.t)^2) G <- cbind(G1,G2,G3,G4) Yhat <- beta.old[1] + (beta.old[2]-beta.old[1]) * (exp.t/(1+exp.t)) e <- Y - Yhat beta.new <- beta.old + solve(t(G) %*% G) %*% t(G) %*% e print(beta.new) diff.beta <- sum((beta.new-beta.old)^2) beta.old <- beta.new } exp.t <- exp(beta.old[3] * (X.t - beta.old[4])) G1 <- 1 - (exp.t/(1+exp.t)) G2 <- exp.t/(1+exp.t) G3 <- (beta.old[2]-beta.old[1]) * (X.t-beta.old[4]) * (exp.t/(1+exp.t)^2) G4 <- (beta.old[2]-beta.old[1]) * (-beta.old[3]) * (exp.t/(1+exp.t)^2) G <- cbind(G1,G2,G3,G4) Yhat <- beta.old[1] + (beta.old[2]-beta.old[1]) * (exp.t/(1+exp.t)) e <- Y - Yhat (SSE <- t(e) %*% e) (MSE <- SSE/(length(Y) - ncol(G))) V.beta <- MSE[1,1] * solve(t(G) %*% G) SE.beta <- sqrt(diag(V.beta)) t.beta <- beta.new/SE.beta P.beta <- 2*(1-pt(abs(t.beta),length(Y)-ncol(G))) print(round(cbind(beta.new,SE.beta,t.beta,P.beta),4)) (CI.beta1 <- beta.new[1] + qt(c(.025,.975),length(Y)-ncol(G))*SE.beta[1]) (CI.beta2 <- beta.new[2] + qt(c(.025,.975),length(Y)-ncol(G))*SE.beta[2]) (CI.beta3 <- beta.new[3] + qt(c(.025,.975),length(Y)-ncol(G))*SE.beta[3]) (CI.beta4 <- beta.new[4] + qt(c(.025,.975),length(Y)-ncol(G))*SE.beta[4]) (h.beta <- 100*(beta.old[2]-beta.old[1]) / beta.old[1]) (H.beta1 <- -100*beta.old[2]/(beta.old[1]^2)) (H.beta2 <- 100/beta.old[1]) H.beta3 <- 0; H.beta4 <- 0 H.beta <- cbind(H.beta1,H.beta2,H.beta3,H.beta4) (SE.h.beta <- sqrt(MSE[1,1] * (H.beta %*% solve(t(G) %*% G) %*% t(H.beta)))) (CI.h.beta <- h.beta + qt(c(.025,.975),length(Y)-ncol(G))*SE.h.beta)