#This is the small model with variation added in at the end. library(arm) library(mvtnorm) #my function to index non-ordered values index.me <- function(v){ x <- sort(unique(v)) index <- rep(0,length(v)) level <- 1:length(x) for(i in 1:length(v)){ for (j in 1:length(x)){ if(v[i]==x[j]) index[i] <- j } } return(index)} #our function to do mle time series estimates ts.mle<-function(data.mat,param.num) { if(!is.matrix(data.mat)){data.mat <- t(as.matrix(data.mat))} out.matrix <- matrix(nrow=dim(data.mat)[1],ncol=param.num+1) for (i in 1:dim(data.mat)[1]){ y <- diff(data.mat[i,]) my.n <- exp(data.mat[i,])-.1 n <- 1:dim(data.mat)[2] - 1 nt <- dim(data.mat)[2] if(param.num==0){ out.matrix[i,1] <- sum(y^2)/ nt c.name <- "sigma" } if (param.num==1){ out.matrix[i,1] <- mean(y) out.matrix[i,2] <- (sum((y-mean(y))^2)) / nt c.name <- c("rhat","sigma") } if(param.num==2){ out.matrix[i,1] <- cov(my.n[n],y)/var(my.n[n]) out.matrix[i,2] <- mean(y) - out.matrix[i,1]*mean(my.n) out.matrix[i,3] <- (sum((y-out.matrix[i,2]-out.matrix[i,1]*my.n[n])^2)) / nt c.name <- c("DD.hat","r.hat","sigma")} } colnames(out.matrix)<- c.name return(out.matrix) } #Here I read in all the data files chiro.log <- read.csv("chirolog.txt") w.means <- read.csv("wmeans.txt") w.cv <- read.csv("wcs.txt") #I need to add lines to strip the read files of row indexes w.means <- as.vector(w.means[,2]) w.cv <- as.vector(w.cv[,2]) chiro.log <- chiro.log[,2:17] chiro.log <- as.matrix(chiro.log) #the function for logistic growth, but can be any function. ft.logis <- function(xx,b){ xx+b[1]+b[2]*exp(xx) } ft1 <- ft.logis nt <- 16 #a metropolis step for sampling a new time series. #Time series length TS <- 16 metnorm <- function(mv1,mv2,vvec){ m1 <- sum(-.5/vvec*(mv1[,1]-mv1[,2])^2) m2 <- sum(-.5/vvec*(mv2[,1]-mv2[,2])^2) a <- exp(m1-m2) z <- runif(1,0,1) xt <- mv2[1,1] if(a>z)xt <- mv1[1,1] return(xt) } #M.means is all the mean water levels. #w.means.i is an index of means to use (i.e. 49 1 - 7) to assign a mean water level value to an estimate #Here is the function to update Beta 1 for water level. b.update.j <- function(y,X,index,my.gamma,J,sigma.b1,sigma.y){ y.temp <- y - X[index,]%*%my.gamma eta.new <- rep(NA,J) for(j in 1:J){ n.j <- sum(index==j) y.bar.j <- mean(y.temp[index==j]) eta.hat.j <- ((n.j/sigma.y^2)*y.bar.j/(n.j/sigma.y^2 + 1/sigma.b1^2)) V.eta.j <- 1/(n.j/sigma.y^2 + 1/sigma.b1^2) eta.new[j] <- rnorm(1,eta.hat.j,sqrt(V.eta.j)) } b1.new.j <- X%*%my.gamma + eta.new return(b1.new.j) } #a version of the function that uses non-nested treatments #This function updates gamma gamma.update <- function(b1,X) { lm.0 <- lm(b1 ~X) g.new <- sim(lm.0,n.sims=1) return(g.new)} sigma.y.update <- function(y,b1,index,n){ sigma.y.new <- sqrt(sum((y-b1[index])^2)/rchisq(1,n-1)) return(sigma.y.new)} sigma.yjk.update <- function(y,y.est){ sigma.y.new <- sqrt(sum((y-y.est)^2)/rchisq(1,length(y)-1)) return(sigma.y.new)} sigma.b1.update <- function(b1,X,my.gamma,J){ sigma.a.new <- sqrt(sum((b1-X%*%my.gamma)^2)/rchisq(1,J-1)) return(sigma.a.new)} #here I'll generate a new time series using the metropolis hastings function I've already written new.mh.timeseries <- function(x,y,ft1,bg,tg,sg) { nt <- length(x) xstar <- rnorm(1,x[1],.05) mv1 <- rbind(c(xstar,y[1]),c(x[2],ft1(xstar,bg))) mv2 <- rbind(c(x[1],y[1]),c(x[2],ft1(x[1],bg))) x[1] <- metnorm(mv1,mv2,c(sg,tg)) #now we're in the middle for(t in 2:(nt-1)){ xstar<-rnorm(1,x[t],.05) mv1 <- rbind(c(xstar,ft1(x[t-1],bg)),c(x[t+1],ft1(xstar,bg))) mv2 <- rbind(c(x[t],ft1(x[t-1],bg)),c(x[t+1],ft1(xstar,bg))) vvec <- c(sg,sg) if(!is.na(y[t])){ mv1 <- rbind(mv1,c(y[t],xstar)) mv2 <- rbind(mv2,c(y[t],x[t])) vvec <- c(vvec,tg) } x[t] <- metnorm(mv1,mv2,vvec) } #last value xstar <- rnorm(1,x[nt],.05) mv1 <- rbind(c(xstar,ft1(x[nt-1],bg)),c(y[nt],xstar)) mv2 <- rbind(c(x[nt],ft1(x[nt-1],bg)),c(y[nt],x[nt])) x[nt] <- metnorm(mv1,mv2,c(sg,tg)) return(x) } #many of these values will be changed to be more generalizable. ts.tau.update <- function(y,x){ v1 <- TS/2 + TS/2 v2 <- .02*((nt/2)-1) + .5*sum((y[!is.na(y)]-x[!is.na(y)])^2) return(1/rgamma(1,v1,v2)) } ts.sig.update<-function(x){ index1 <- 1:TS-1 index2 <- 2:TS w1 <- TS/10 + TS/2 w2 <- .1*((TS/(TS*.5))-1)+ .5*sum((x[index2]-ft1(x[index1],bg))^2) return(1/rgamma(1,w1,w2)) } #here is the function to draw a new pair of beta's from the time series. draw.betas <- function(x){ nt <- length(x) index1 <- 1:nt-1 prior.IVB <- diag(.00001,2) prior.B <- c(.01,.01) xmat <- cbind(rep(1,(nt-1)),exp(x[index1])) difx <- as.matrix(diff(x)) bigv <- solve(1/sg*crossprod(xmat)+prior.IVB) smallv <- 1/sg*crossprod(xmat,difx)+t(prior.B%*%prior.IVB) bg <- rmvnorm(1,bigv%*%smallv,bigv) return(bg) } #the normal sample I actually want to use norm.samp <- function(y){ y.temp <- y n <- length(y) ybar <- mean(y.temp) Ss <- sum((y.temp-ybar)^2) sig <- Ss/rchisq(1,length(y)-1) mu <- rnorm(1,ybar,sqrt(sig)/sqrt(n)) return(list(mu=mu,sig=sig))} #I'm not sure what this norm samp is #norm.samp <- function(y,X.f1,index.f1,my.gamma.f1,X.f2,index.f2,my.gamma.f2){ #y.temp <- y - X.f1[index.f1,]%*%my.gamma.f1 - X.f2[index.f2,]%*%my.gamma.f2 #ybar <- mean(y.temp) #Ss <- sum((y.temp-ybar)^2) #sig <- Ss/rchisq(1,length(y)-1) #mu <- rnorm(1,ybar,sqrt(sig)/sqrt(n)) #return(list(mu=mu,sig=sig))} #X is a design matrix of the different treatment levels with a constant in the first column #ok, now to put it all together for a little test X <- cbind(w.cv,index.me(w.cv),w.means,index.me(w.means)) X.f1 <- cbind(rep(1,7),sort(unique(X[,1]))-mean(sort(unique(X[,1])))) X.f1.labels <- sort(unique(X[,1])) index.f1 <- X[,2] X.f2 <- cbind(rep(1,7),sort(unique(X[,3]))-mean(sort(unique(X[,3])))) index.f2 <- X[,4] X.f2.labels <- sort(unique(X[,3])) J <- 7 #the number of groups K <- 7 N<- J*K TS <- 16 test.mle <- ts.mle(chiro.log,2) #save.image(file="Chiro.Analysis-3mods.rdata") ngibbs <- 10000 sg <- runif(1) tg <- runif(1) xgibbs <- array(dim=c(ngibbs,TS,N)) bgibbs.est <- array(dim=c(ngibbs,2,N)) tgibbs <- array(dim=c(ngibbs,1,N)) sgibbs <- array(dim=c(ngibbs,1,N)) b0.gibbs.f1 <- matrix(nrow=ngibbs,ncol=J) b1.gibbs.f1 <- matrix(nrow=ngibbs,ncol=J) sig.b1 <- rep(NA,ngibbs) sig.b0 <- rep(NA,ngibbs) sig.k <- rep(NA,ngibbs) sig.b1.group.f1 <- rep(NA,ngibbs) sig.b0.group.f1 <- rep(NA,ngibbs) sig.b1.indiv.f2 <- rep(NA,ngibbs) sig.b0.indiv.f2 <- rep(NA,ngibbs) sig.b1.indiv.f1 <- rep(NA,ngibbs) sig.b0.indiv.f1 <- rep(NA,ngibbs) b0.gibbs.f1 <- matrix(nrow=ngibbs,ncol=J) b1.gibbs.f1 <- matrix(nrow=ngibbs,ncol=J) sig.b1.group.f1 <- rep(NA,ngibbs) sig.b1 <- rep(NA,ngibbs) sig.b0 <- rep(NA,ngibbs) sig.k <- rep(NA,ngibbs) sig.b0.group.f1 <- rep(NA,ngibbs) mu.b0 <- rep(NA,ngibbs) mu.b1 <- rep(NA,ngibbs) sig.b1.indiv.f1 <- rep(NA,ngibbs) sig.b0.indiv.f1 <-rep(NA,ngibbs) sig.b1.indiv.f2 <- rep(NA,ngibbs) sig.b0.indiv.f2 <-rep(NA,ngibbs) #Since we don't have any group level parameters at the moment, we'll just have them #be random effects #b0.gamma.f1 <- matrix(nrow=ngibbs,ncol=2) #b1.gamma.f1 <- matrix(nrow=ngibbs,ncol=2) b0.gibbs.f2 <- matrix(nrow=ngibbs,ncol=K) b1.gibbs.f2 <- matrix(nrow=ngibbs,ncol=K) sig.b1.group.f2 <- rep(NA,ngibbs) sig.b0.group.f2 <-rep(NA,ngibbs) #b0.gamma.f2 <- matrix(nrow=ngibbs,ncol=2) #b1.gamma.f2 <- matrix(nrow=ngibbs,ncol=2) b1.gibbs.adj.f1<-matrix(nrow=ngibbs,ncol=J) b0.gibbs.adj.f1<-matrix(nrow=ngibbs,ncol=J) b1.gibbs.adj.f2<-matrix(nrow=ngibbs,ncol=K) b0.gibbs.adj.f2<-matrix(nrow=ngibbs,ncol=K) b0.gamma.f1 <- matrix(nrow=ngibbs,ncol=2) b1.gamma.f1 <- matrix(nrow=ngibbs,ncol=2) b0.gibbs.f2 <- matrix(nrow=ngibbs,ncol=K) b1.gibbs.f2 <- matrix(nrow=ngibbs,ncol=K) sig.b1.group.f2 <- rep(NA,ngibbs) sig.b0.group.f2 <-rep(NA,ngibbs) b0.gamma.f2 <- matrix(nrow=ngibbs,ncol=2) b1.gamma.f2 <- matrix(nrow=ngibbs,ncol=2) b1.gibbs.adj.f1<-matrix(nrow=ngibbs,ncol=J) b0.gibbs.adj.f1<-matrix(nrow=ngibbs,ncol=J) b1.gibbs.adj.f2<-matrix(nrow=ngibbs,ncol=K) b0.gibbs.adj.f2<-matrix(nrow=ngibbs,ncol=K) #We have to put in initial parameters in each data array #this inputs the actual data in the initial time series for(i in 1:N) { xgibbs[1,,i]<- chiro.log[i,] bgibbs.est[1,,i] <- test.mle[i,c(2,1)] sgibbs[1,1,i] <- test.mle[i,3] tgibbs[1,1,i] <- .2 } sig.b1.group.f1[1] <- runif(1,0,1) sig.b1.indiv.f1[1] <- runif(1,0,1) sig.b0.group.f1[1] <- runif(1,0,1) sig.b0.indiv.f1[1] <- runif(1,0,1) sig.b1.group.f2[1] <- runif(1,0,1) sig.b1.indiv.f2[1] <- runif(1,0,1) sig.b0.group.f2[1] <- runif(1,0,1) sig.b0.indiv.f2[1] <- runif(1,0,1) b0.gibbs.f1[1,] <- rnorm(J,tapply(test.mle[,2],index.f1,mean),.01) b1.gibbs.f1[1,] <- -1*rexp(J,-1/tapply(test.mle[,1],index.f1,mean)) b1.gamma.f1[1,] <- c(-.0004,-.0004) b0.gamma.f1[1,] <- c(1,1) b1.gibbs.f2[1,] <- rnorm(K,tapply(test.mle[,2],index.f1,mean),.01) b0.gibbs.f2[1,] <- -1*rexp(K,-1/tapply(test.mle[,1],index.f1,mean)) b1.gibbs.adj.f1[1,] <- b1.gibbs.f1[1,] - mean(b1.gibbs.f1[1,]) b0.gibbs.adj.f1[1,] <- b0.gibbs.f1[1,] - mean(b0.gibbs.f1[1,]) b1.gibbs.adj.f2[1,] <- b1.gibbs.f2[1,] - mean(b1.gibbs.f2[1,]) b0.gibbs.adj.f2[1,] <- b0.gibbs.f2[1,] - mean(b0.gibbs.f2[1,]) sig.b1.group.f2[1] <- .0000003 sig.b1[1] <- .00003 sig.b0.group.f2[1] <- .00000003 sig.b0[1] <- .000003 mu.b0[1] <- rnorm(1,mean(test.mle[,2]),.01) mu.b1[1] <- mean(test.mle[,1]) b1.gibbs.f2[1,] <- tapply(test.mle[,1],index.f2,mean) b0.gibbs.f2[1,] <- tapply(test.mle[,2],index.f2,mean) b1.gamma.f2[1,] <- c(-.0004,-.0004) b0.gamma.f2[1,] <- c(.0001,.0001) act.est.b0 <- matrix(nrow=ngibbs,ncol=N) act.est.b1 <- matrix(nrow=ngibbs,ncol=N) alt.est.b1 <- matrix(nrow=ngibbs,ncol=N) alt.est.b0 <- matrix(nrow=ngibbs,ncol=N) act.est.b1[1,] <-test.mle[,1] act.est.b0[1,] <- test.mle[,2] alt.est.b1[1,] <-test.mle[,1] alt.est.b0[1,] <- test.mle[,2] #Here I will set up the array that will hold my data full data structure which can then be treated as a "bugs" object for (i in 1:(ngibbs-1)){ for(z in 1:N){ #new.mh.timeseries <- function(x,y,ft1,bg,tg,sg) #draw betas to use to estimate; # if(is.nan(rexp(1,1/(-1*(mu.b1[i]+b1.gibbs.adj.f2[i,index.f2[z]]+b1.gibbs.adj.f1[i,index.f1[z]]))))==1){ # my.b1 <- 0} else {my.b1 <- -1*rexp(1,1/(-1*(mu.b1[i]+b1.gibbs.adj.f2[i,index.f2[z]]+b1.gibbs.adj.f1[i,index.f1[z]]))) } my.b1 <- mu.b1[i]+b1.gibbs.adj.f2[i,index.f2[z]]+b1.gibbs.adj.f1[i,index.f1[z]] mh.beta <- c((mu.b0[i]+b0.gibbs.adj.f2[i,index.f2[z]]+b0.gibbs.adj.f1[i,index.f1[z]]),my.b1) act.est.b0[i+1,z] <- mh.beta[1] act.est.b1[i+1,z] <- mh.beta[2] x <- new.mh.timeseries(xgibbs[i,,z],xgibbs[1,,z],ft1,mh.beta,tgibbs[i,1,z],sgibbs[i,1,z]) bg <- draw.betas(x) xgibbs[i+1,,z] <- x bgibbs.est[i+1,,z] <- bg sg <- ts.sig.update(xgibbs[i,,z]) #NOTE THIS CHANGE FOR THE MISSING VALUES OF THE CHIRO SERIES! #SWITCH BACK TO xgibbs[1,,z] tg <- ts.tau.update(chiro.log[z,],xgibbs[i,,z]) tgibbs[i+1,1,z] <- tg sgibbs[i+1,1,z] <- sg } b0.data <- bgibbs.est[i,1,] b1.data <- bgibbs.est[i,2,] b1.gibbs.f1[i+1,] <- b.update.j(b1.data,X.f1,index.f1,b1.gamma.f1[i,],J,sig.b1.group.f1[i],sig.b1.indiv.f1[i]) b0.gibbs.f1[i+1,] <- b.update.j(b0.data,X.f1,index.f1,b0.gamma.f1[i,],J,sig.b0.group.f1[i],sig.b0.indiv.f1[i]) b1.gamma.f1[i+1,] <- gamma.update(b1.gibbs.adj.f1[i,],X.f1)$coef b0.gamma.f1[i+1,] <- gamma.update(b0.gibbs.adj.f1[i,],X.f1)$coef sig.b1.group.f1[i+1] <- sigma.b1.update(b1.gibbs.f1[i,],X.f1,b1.gamma.f1[i,],J) sig.b1.indiv.f1[i+1] <- sigma.y.update(b1.data,b1.gibbs.f1[i,],index.f1,49) sig.b0.group.f1[i+1] <- sigma.b1.update(b0.gibbs.f1[i,],X.f1,b0.gamma.f1[i,],J) sig.b0.indiv.f1[i+1] <- sigma.y.update(b0.data,b0.gibbs.f1[i,],index.f1,49) # mu.b1.t <- rexp(1,1/mean(-1*b1.data)) # mu.b1[i+1] <- -1*mu.b1.t mu.b1[i+1] <- norm.samp(b1.data)$mu sig.b1[i+1] <- sqrt(norm.samp(b1.data)$sig) mu.b0[i+1] <- norm.samp(b0.data)$mu sig.b0[i+1] <- sqrt(norm.samp(b0.data)$sig) b1.gibbs.f2[i+1,] <- b.update.j(b1.data,X.f2,index.f2,b1.gamma.f2[i,],K,sig.b1.group.f2[i],sig.b1.indiv.f2[i]) b0.gibbs.f2[i+1,] <- b.update.j(b0.data,X.f2,index.f2,b0.gamma.f2[i,],K,sig.b0.group.f2[i],sig.b0.indiv.f2[i]) b1.gamma.f2[i+1,] <- gamma.update(b1.gibbs.adj.f2[i,],X.f2)$coef b0.gamma.f2[i+1,] <- gamma.update(b0.gibbs.adj.f2[i,],X.f2)$coef sig.b1.group.f2[i+1] <- sigma.b1.update(b1.gibbs.f2[i,],X.f2,b1.gamma.f2[i,],K) sig.b1.indiv.f2[i+1] <- sigma.y.update(b1.data,b1.gibbs.f2[i,],index.f2,49) sig.b0.group.f2[i+1] <- sigma.b1.update(b0.gibbs.f2[i,],X.f2,b0.gamma.f2[i,],K) sig.b0.indiv.f2[i+1] <- sigma.y.update(b0.data,b0.gibbs.f2[i,],index.f2,49) b1.gibbs.adj.f1[i+1,] <- b1.gibbs.f1[i,] - mean(b1.gibbs.f1[i,]) b0.gibbs.adj.f1[i+1,] <- b0.gibbs.f1[i,] - mean(b0.gibbs.f1[i,]) b1.gibbs.adj.f2[i+1,] <- b1.gibbs.f2[i,] - mean(b1.gibbs.f2[i,]) b0.gibbs.adj.f2[i+1,] <- b0.gibbs.f2[i,] - mean(b0.gibbs.f2[i,]) cat(i, " ") } X.f2.labels <- sort(unique(X[,3])) X.f1.labels <- sort(unique(X[,1])) c.names.b1 <- c(paste("b1.",X.f1.labels,".CV"),"sig.b1.CV","gamma.f1.intercept","gamma.f1.slope",paste("b1.",round(X.f2.labels,0),".WL"),"sig.b1.WL","gamma.b1.intercept","gamma.b1.slope") c.names.b0 <- c(paste("b0.",X.f1.labels,".CV"),"sig.b0.CV","gamma.f1.intercept","gamma.f1.slope",paste("b0.",round(X.f2.labels,0),".WL"),"sig.b0.WL","gamma.f2.intercept","gamma.f2.slope") c.names.est.b1 <- paste("b1",X.f1.labels[index.f1],round(X.f2.labels[index.f2],0)) c.names.est.b0 <- paste("b0",X.f1.labels[index.f1],round(X.f2.labels[index.f2],0)) b1.output <- matrix(nrow=ngibbs,ncol=K+J+6) b0.output <- matrix(nrow=ngibbs,ncol=K+J+6) for(i in 1:ngibbs){ b1.output[i,] <- c(b1.gibbs.f1[i,],sig.b1.group.f1[i],b1.gamma.f1[i,],b1.gibbs.f2[i,],sig.b1.group.f2[i],b1.gamma.f2[i,]) b0.output[i,] <- c(b0.gibbs.f1[i,],sig.b0.group.f1[i],b0.gamma.f1[i,],b0.gibbs.f2[i,],sig.b0.group.f2[i],b0.gamma.f2[i,]) } colnames(act.est.b1) <- c.names.est.b1 colnames(act.est.b0) <- c.names.est.b0 colnames(b1.output) <- c.names.b1 colnames(b0.output) <- c.names.b0 chiro.b1.output <- b1.output chiro.b0.output <- b0.output chiro.act.est.b1 <- act.est.b1 chiro.act.est.b0 <- act.est.b0 chiro.xgibbs <- xgibbs chiro.bgibbs <- bgibbs.est chiro.sgibbs <- sgibbs chiro.tgibbs <- tgibbs chiro.sig.b0 <- sig.b0 chiro.sig.b1 <- sig.b1 chiro.mu.b1 <- mu.b1 chiro.mu.b0 <- mu.b0 save.image(file="Chiro.Analysis.rdata")