options(object.size=10^100) ###FIRST DEFINE FUNCTION ###doit() fits using local likelihood, window makes the weights, ###geninv returns the general inverse of a matrix ##THIS IS FUNCTION THAT FITS THE LOCAL-LIKELIHOOD ##The program creates a plot and calls it plot04.ps doit <- function(){ res <- rep(0,M) for(t0 in o){ print(t0) if(t0 < (ws+1)/2) { tt <- 1:ws ww <- sapply(c((t0-1):0,1:(ws-t0))/(ws-t0+1),function(x) (1-x^3)^3) middle <- t0 } else{ if(t0 > (M - (ws+1)/2 + 1)){ tt <- (M-ws+1):M ww <- sapply(c((t0-(M-ws+1)):1,0:(M-t0))/(t0-(M-ws)), function(x) (1-x^3)^3) middle <- ws - M + t0 } else{ tt <- (t0 - (ws-1)/2):(t0 + (ws-1)/2) ww <- sapply(c(((ws-1)/2):0,1:((ws-1)/2))/((ws+1)/2), function(x) (1-x^3)^3) middle <- (ws+1)/2 } } xx <- rep(Age[tt], N[tt]) yy <- c() for(h in tt) yy <- c(yy, rep(1, Y[h]), rep(0, N[h] - Y[h])) ww <- rep(ww, N[tt]) assign("xx",xx,frame=1) X <- sapply(0:degree, function(k) xx^k) fit1 <- gam(yy ~ X, family = "binomial", control = glm.control(epsilon = 1e-008), weights = ww) res[t0] <- fit1\$fitted[(cumsum(N[tt]))[middle]] } lines(Age,res,lty=LTY) } det <- function(M) prod(svd(M, nu = 0, nv = 0)\$d) geninv <- function(M){ aux <- svd(M) auxd <- aux\$d auxd[auxd!=0] <- 1/auxd aux\$v%*%diag(auxd,nrow=length(auxd))%*%t(aux\$u) } Window <- function(n, h=1) { aux <- seq(-0.999, 0.999, len = round(n * h)) Naux1 <- round(n - length(aux))/2 Naux2 <- n - length(aux) - Naux1 c(rep(0, Naux1), (1 - abs(aux)^3)^3, rep(0, Naux2)) } Data <- matrix(scan("trevor.dat"),byrow=T,ncol=4) ###entries are survive,age,year of operation (as in 19??) and Number of nodes. ###See Tibs and Hast 1987 Age <- Data[,2] #N <- Data[,2] Y <- tapply(Data[,1],Age,sum) N <- tapply(Data[,2],Age,length) Age <- unique(Age) N0 <- 25 NS <- seq(4,24) Ks <- c(0:2) Criteria <- array(0,dim=c(6,length(NS),length(Ks))) for(i in 1:length(NS)){ if(NS[i] <= 24) {o <- (N0 - NS[i]):(N0 + NS[i]);M <- length(o)} else {o <- 1:49;M <- 2 * NS[i] + 1} xx <- rep(Age[o],N[o]) yy <- c(); for(h in o) yy <- c(yy,rep(1,Y[h]),rep(0,N[h]-Y[h])) ww <- Window(M,1); if(length(ww)>length(o)) ww <- ww[((M-1)/2 - 24):((M-1)/2 + 24)] ww <- rep(ww,N[o]) for(j in 1:length(Ks)){ X <- sapply(0:Ks[j],function(k) scale01(xx,-1,1)^k) fit1 <- glm(yy~X-1,family="binomial",weights=ww) p <- fit1\$fitted AUX <- sweep(X,1,ww,FUN="*") AUX <- sweep(AUX,1,p*(1-p),FUN="*") AUX0 <- t(XX)%*%XX AUX1 <- t(AUX)%*%X AUX <- sweep(AUX,1,ww,FUN="*") AUX2 <- t(AUX)%*%X NP <- sum(diag(geninv(AUX1) %*% AUX2)) l0 <- sum(ww*(yy*log(p) + (1-yy)*log(1-p))) Criteria[1,i,j] <- (-2*l0 + 2*j)/length(N) Criteria[2,i,j] <- (-2*l0 + 2*NP)/sum(ww) Criteria[3,i,j] <- (-2*l0 + j*log(det(AUX0)))/length(N) Criteria[4,i,j] <- (-2*l0 + log(det(AUX1)))/sum(ww) Criteria[5,i,j] <- (-2*l0 + j*log(det(AUX0)) + 2*j)/length(N) Criteria[6,i,j] <- (-2*l0 + log(det(AUX1)) + 2*NP)/sum(ww) } } postscript("Plots/plot04.ps") Title <- c("WAIC","WBIC","WCAICF") par(mfrow=c(2,2), mar=c(3.5,2,2,.2),mgp=c(2,.5,0)) for(i in 1:3){ plot(2*NS+1,Criteria[i*2,,1],ylim=range(Criteria[i*2,,]),type="l",main=Title[i],yaxt="n",xlab="Window Size in Years") lines(2*NS+1,Criteria[i*2,,2],lty=2) lines(2*NS+1,Criteria[i*2,,3],lty=3) if(i==1) legend(22,1.27,legend=c("constant","line","parabola"),lty=c(1,2,3),cex=.8) } aux <- matrix(0,6,2) for(k in 1:6){ x <- order(as.vector(Criteria[k,,]))[1] i <- ceiling(x/length(NS)) aux[k,]<-c(i,NS[x - (i-1)*length(NS)]) } M <- length(N) o <- 1:M plot(Age,Y/N,xlab="Age in Years",ylab="Percent Survived") for(i in c(2,4,6)){ degree <- aux[i,1]-1;ws <- 2*aux[i,2]+1;LTY <- i/2 doit() } legend(32,.32,legend=c("WAIC","WBIC","WCAICF"),lty=c(1,2,3),cex=.8) dev.off()