= Bits of SPlus Code for analyses in 2003 =

=== get.roi.data.linear ===

{{{
get.roi.data.linear <- function(subj, cond, colnames)
{
	file <- paste("C:\\Documents and Settings\\user\\My Documents\\People\\Matthew Brett\\Time Series\\roi_data_linear_",
		subj, "_", cond, ".wk1", sep = "")
	frame <- paste("data.", subj, ".", cond, sep = "")
	print(file)
	import.data(FileName = file, FileType = "Lotus", ColNames = colnames, 
		TargetStartCol = "1", DataFrame = frame, StartCol = "1", EndCol = "END", 
		StartRow = "1", EndRow = "END", Delimiters = ", \\t", SeparateDelimiters = 
		F, PageNumber = "1", RowNameCol = "", StringsAsFactors = "Auto", 
		VLabelAsNumber = F, DatesAsDoubles = F, Filter = "")
	source(text = paste("names(", frame, ")<-colnames", sep = ""))
	guiClose("DataFrame", frame)
}
}}}

=== pacf ===

{{{
pacf <- function(s) {
	acf(s, type = "partial")
}
}}}

=== loopMLlinear ===

{{{
loopMLlinear <- function() {
	graphsheet(page=T)
	results <- c()
	for (subj in subjects) { 
	for (cond in conditions) {
		get.roi.data.linear(subj,cond,colnames)
#		dev.ask()
		dataname <- paste("data.", subj, ".", cond, sep = "")
		source(text=paste("attach(", dataname,")",sep=""))
		result <- actionML(dataname)
		results <- c(results, result)
		source(text=paste("detach(\"", dataname,"\")",sep=""))
	}
}
return(results=results) 
}
}}}



=== innovations ===

{{{
innovations <- function(x,theta=c()) {
	n <- length(x)
	p <- length(theta)
	common.times <- (p+1):n
	y <- x[common.times]
	if (p==0) { return(y) }
	for (k in 1:p) { 	y <- y - theta[k]*x[common.times-k] }
	return(y)
}
}}}

=== meld ===

{{{
meld <- function(su,co) {
	paste("data.",su,".",co,sep="")
}
}}}

=== gather_phiML ===

{{{
phisML <- c()
for (subj in 1:12) {
	for (cond in 1:7) {
		datastart <- ((subj-1)*7+(cond-1))*6
		phisML <- c(phisML, coef(resultsML[[datastart+4]]$modelStruct$corStruct,F))
	}
}

phisMLt <- t(matrix(phisML,nrow=3))

export.data("phisMLt","phisMLlinear.xls","EXCEL")
}}}

=== gather_phis ===

{{{
# coef(results[[10]]$modelStruct$corStruct,F)

phis <- c()
for (subj in 1:12) {
	for (cond in 1:7) {
		datastart <- ((subj-1)*7+(cond-1))*12
		phis <- c(phis, coef(results[[datastart+10]]$modelStruct$corStruct,F))
	}
}
}}}

=== gather_coeffsML ===

{{{
coeffsML <- c()
for (subj in 1:12) {
	for (cond in 1:7) {
		datastart <- ((subj-1)*7+(cond-1))*6
		coeffsML <- c(coeffsML, resultsML[[datastart+1]]$coefficients[2:3,1:2])
		coeffsML <- c(coeffsML, resultsML[[datastart+4]]$tTable[2:3,1:2])
	}
}
}}}


=== gather_coeffs ===

{{{
coeffs <- c()
for (subj in 1:12) {
	for (cond in 1:7) {
		datastart <- ((subj-1)*7+(cond-1))*6
		coeffs <- c(coeffs, results[[datastart+1]]$coefficients[2:3,1:2])
		coeffs <- c(coeffs, results[[datastart+4]]$tTable[2:3,1:2])
	}
}

coeffs <- t(matrix(coeffs,nrow=8))

export.data("coeffs","coeffsOLS+MLlinear.xls","EXCEL")
coeffsML <- t(matrix(coeffsML,nrow=2))

export.data("coeffsML","coeffsMLlinear.xls","EXCEL")
}}}


=== gather_sigmaML ===

{{{
sigmaML <- c()
for (subj in 1:12) {
	for (cond in 1:7) {
		datastart <- ((subj-1)*7+(cond-1))*6
		sigmaML <- c(sigmaML, resultsML[[datastart+1]]$sigma)
		sigmaML <- c(sigmaML, resultsML[[datastart+4]]$sigma)
	}
}

sigmaML <- t(matrix(sigmaML,nrow=2))

export.data("sigmaML","sigmaMLlinear.xls","EXCEL")
}}}


=== print.summary.gls ===

{{{
function(x, verbose = FALSE, digits = .Options$digits, ...)
{
	dd <- x$dims
	verbose <- verbose || attr(x, "verbose")
	mCall <- x$call
	if(inherits(x, "gnls")) {
		cat("Generalized nonlinear least squares fit\n")
	}
	else {
		cat("Generalized least squares fit by ")
		cat(ifelse(x$method == "REML", "REML\n", "maximum likelihood\n"))
	}
	cat("  Model:", deparse(as.vector(mCall$model)), "\n")
	cat("  Data:", deparse(mCall$data), "\n")
	if(!is.null(mCall$subset)) {
		cat("  Subset:", deparse(asOneSidedFormula(mCall$subset)[[2]]), "\n")
	}
	print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, row.names = " "))
	if(verbose) {
		cat("Convergence at iteration:", x$numIter, "\n")
	}
	if(length(x$modelStruct)) {
		cat("\n")
		print(summary(x$modelStruct))
	}
	cat("\nCoefficients:\n")
	xtTab <- as.data.frame(x$tTable)
	wchPval <- match("p-value", names(xtTab))
	for(i in names(xtTab)[ - wchPval]) {
		xtTab[, i] <- format(zapsmall(xtTab[, i]))
	}
	xtTab[, wchPval] <- format(round(xtTab[, wchPval], 4))
	if(any(wchLv <- (as.double(levels(xtTab[, wchPval])) == 0))) {
		levels(xtTab[, wchPval])[wchLv] <- "<.0001"
	}
	row.names(xtTab) <- dimnames(x$tTable)[[1]]
	print(xtTab)
	if(nrow(x$tTable) > 1) {
		corr <- x$corBeta
		class(corr) <- "correlation"
		print(corr, title = "\n Correlation:", ...)
	}
	cat("\nStandardized residuals:\n")
	print(x$residuals)
	cat("\n")
	cat("Residual standard error:", format(x$sigma), "\n")
	cat("Degrees of freedom:", dd[["N"]], "total;", dd[["N"]] - dd[["p"]], "residual\n")
}
}}}

=== dwstat ===

{{{
dwstat <- function(resids, S, X, mu, max.lag = 1, simulate = TRUE, reps = 1000, method = c(
	"resample", "normal"), alternative = c("two.sided", "positive", 
	"negative"))
{
	method <- match.arg(method)
	alternative <- if(max.lag == 1) match.arg(alternative) else 
			"two.sided"
#	resids <- residuals(model)
	if(any(is.na(resids)))
		stop("residuals include missing values")
	n <- length(resids)
	r <- dw <- rep(0, max.lag)
	den <- sum(resids^2)
	for(lag in 1:max.lag) {
		dw[lag] <- (sum((resids[(lag + 1):n] - resids[1:(n - 
			lag)])^2))/den
		r[lag] <- (sum(resids[(lag + 1):n] * resids[1:(n - lag)
			]))/den
	}
	if(!simulate) {
		result <- list(r = r, dw = dw)
		class(result) <- "durbin.watson"
		result
	}
	else {
#		S <- summary(model)$sigma
#		X <- model.matrix(model)
#		mu <- fitted.values(model)
		Y <- if(method == "resample") matrix(sample(resids, n * 
				reps, replace = TRUE), n, reps) + matrix(mu, 
				n, reps) else matrix(rnorm(n * reps, 0, S), n,
				reps) + matrix(mu, n, reps)
		E <- residuals(lm(Y ~ X - 1))
		DW <- apply(E, 2, durbin.watson, max.lag = max.lag)
		if(max.lag == 1)
			DW <- rbind(DW)
		p <- rep(0, max.lag)
		if(alternative == "two.sided") {
			for(lag in 1:max.lag) {
				p[lag] <- (sum(dw[lag] < DW[lag,  ]))/reps
				p[lag] <- 2 * (min(p[lag], 1 - p[lag]))
			}
		}
		else if(alternative == "positive") {
			for(lag in 1:max.lag) {
				p[lag] <- (sum(dw[lag] > DW[lag,  ]))/reps
			}
		}
		else {
			for(lag in 1:max.lag) {
				p[lag] <- (sum(dw[lag] < DW[lag,  ]))/reps
			}
		}
		result <- list(r = r, dw = dw, p = p, alternative = 
			alternative)
		class(result) <- "durbin.watson"
		result
	}
}
}}}

=== durbin.watson.lm ===

{{{
durbin.watson.lm <- function(model, max.lag = 1, simulate = TRUE, reps = 1000, method = c(
	"resample", "normal"), alternative = c("two.sided", "positive", 
	"negative"))
{
	method <- match.arg(method)
	alternative <- if(max.lag == 1) match.arg(alternative) else 
			"two.sided"
	resids <- residuals(model)
	if(any(is.na(resids)))
		stop("residuals include missing values")
	n <- length(resids)
	r <- dw <- rep(0, max.lag)
	den <- sum(resids^2)
	for(lag in 1:max.lag) {
		dw[lag] <- (sum((resids[(lag + 1):n] - resids[1:(n - 
			lag)])^2))/den
		r[lag] <- (sum(resids[(lag + 1):n] * resids[1:(n - lag)
			]))/den
	}
	if(!simulate) {
		result <- list(r = r, dw = dw)
		class(result) <- "durbin.watson"
		result
	}
	else {
		S <- summary(model)$sigma
		X <- model.matrix(model)
		mu <- fitted.values(model)
		Y <- if(method == "resample") matrix(sample(resids, n * 
				reps, replace = TRUE), n, reps) + matrix(mu, 
				n, reps) else matrix(rnorm(n * reps, 0, S), n,
				reps) + matrix(mu, n, reps)
		E <- residuals(lm(Y ~ X - 1))
		DW <- apply(E, 2, durbin.watson, max.lag = max.lag)
		if(max.lag == 1)
			DW <- rbind(DW)
		p <- rep(0, max.lag)
		if(alternative == "two.sided") {
			for(lag in 1:max.lag) {
				p[lag] <- (sum(dw[lag] < DW[lag,  ]))/reps
				p[lag] <- 2 * (min(p[lag], 1 - p[lag]))
			}
		}
		else if(alternative == "positive") {
			for(lag in 1:max.lag) {
				p[lag] <- (sum(dw[lag] > DW[lag,  ]))/reps
			}
		}
		else {
			for(lag in 1:max.lag) {
				p[lag] <- (sum(dw[lag] < DW[lag,  ]))/reps
			}
		}
		result <- list(r = r, dw = dw, p = p, alternative = 
			alternative)
		class(result) <- "durbin.watson"
		result
	}
}
}}}

=== botch ===

{{{
print.summary.corStruct <-
  function(x, ...)
{
  class(x) <- attr(x, "oClass")
  stNam <- attr(x, "structName")
  form <- formula(x)
  if (!is.null(stNam)) {
    cat(paste("Correlation Structure: ", stNam, "\n", sep = ""))
  }
  if (!is.null(form)) {
    cat(paste(" Formula:", deparse(as.vector(form)),"\n"))
  }
  cat(" Parameter estimate(s):\n")
  print(coef(x, FALSE))
#  return(coef(x, FALSE))
}
}}}

=== ar_4steps ===

{{{
par(mfrow=c(4,1))
d10.ols <- gls(bold~hrf,roi.data.990151.10,method="ML")
summary(d10.ols)
pacf(d10.ols$residuals)
d10.AR1 <- gls(bold~hrf,roi.data.990151.10,correlation=corARMA(p=1),method="ML")
summary(d10.AR1)
pacf(d10.AR1$residuals)
d10.AR2 <- gls(bold~hrf,roi.data.990151.10,correlation=corARMA(p=2),method="ML")
summary(d10.AR2)
pacf(d10.AR2$residuals)
d10.AR3 <- gls(bold~hrf,roi.data.990151.10,correlation=corARMA(p=3),method="ML")
summary(d10.AR3)
pacf(d10.AR3$residuals)
cor(d10.AR3$residuals,lag(d10.AR3$residuals,1))
}}}


=== actionML ===

{{{
actionML <- function(dataname) {
	par(mfrow = c(5, 3))
	tsplot(bold, type = "l")
	title(paste(dataname,": Raw BOLD",sep=""))

# First for the OLS

	ols <- lm(bold~hrf+dhrf)
	errors0 <- residuals(ols)
	fitted0 <- fitted.values(ols)
	tsplot(fitted0)
	title(paste(dataname,": Fitted(OLS)",sep=""))
	tsplot(errors0)
	title(paste(dataname,": Residuals(OLS)",sep=""))
		auto0 <- coef(ols)
	output0 <- acf(errors0, type = "partial", plot = T)
	dw0 <- dwstat(errors0, summary(ols)$sigma, model.matrix(ols),
				fitted.values(ols),max.lag=8)
	plot(cbind(dw0$r,dw0$p),main= 'Durbin Watson for OLS', xlab= 'Autocorrelation', ylab= 'P value')
	
# Now try AR3	
	
	ar3 <- gls(bold~hrf+dhrf,correlation=corARMA(p=3),method="ML")
	auto3 <- coef(ar3$modelStruct$corStruct, F)
	errors3 <- innovations(ar3$resid,auto3)
	current <- as.data.frame(errors3)
	attach(current)
	tsplot(errors3)
	title(paste(dataname,": Innovations(AR3)",sep=""))
	output3 <- acf(errors3, type = "partial", plot = T)
	lm3 <- lm(errors3~1)
	errors3 <- residuals(lm3)
	sigma3 <- summary(lm3)$sigma
	modelmatrix3 <- eval(model.matrix(lm3))
	fittedvalues3 <- fitted.values(lm3)
	dw3 <- dwstat(errors3, sigma3, modelmatrix3, fittedvalues3, max.lag=8)
	plot(cbind(dw3$r,dw3$p),main= 'Durbin Watson for AR3', xlab= 'Autocorrelation', ylab= 'P value')
	detach("current")
	
	return(
		AR0=summary(ols), dw0=dw0, Cols=coefficients(ols),
		AR3=summary(ar3), dw3=dw3, Car3=coefficients(ar3)
	)
}
}}}


=== action ===

{{{
action <- function(dataname) {
	par(mfrow = c(5, 3))
	tsplot(bold, type = "l")
	title(paste(dataname,": Raw BOLD",sep=""))

# First for the OLS

	ols <- lm(bold~hrf+dhrf)
	errors0 <- residuals(ols)
	fitted0 <- fitted.values(ols)
	tsplot(fitted0)
	title(paste(dataname,": Fitted(OLS)",sep=""))
	tsplot(errors0)
	title(paste(dataname,": Residuals(OLS)",sep=""))
		auto0 <- coef(ols)
	output0 <- acf(errors0, type = "partial", plot = T)
	dw0 <- dwstat(errors0, summary(ols)$sigma, model.matrix(ols),
				fitted.values(ols),max.lag=8)
	plot(cbind(dw0$r,dw0$p),main= 'Durbin Watson for OLS', xlab= 'Autocorrelation', ylab= 'P value')
	
# Now try AR1	
	
	ar1 <- gls(bold~hrf+dhrf,correlation=corARMA(p=1))
	auto1 <- coef(ar1$modelStruct$corStruct, F)
	errors1 <- innovations(ar1$resid,auto1)
	current <- as.data.frame(errors1)
	attach(current)
	tsplot(errors1)
	title(paste(dataname,": Innovations(AR1)",sep=""))
	output1 <- acf(errors1, type = "partial", plot = T)
	lm1 <- lm(errors1~1)
	errors1 <- residuals(lm1)
	sigma1 <- summary(lm1)$sigma
	modelmatrix1 <- model.matrix(lm1)
	fittedvalues1 <- fitted.values(lm1)
	dw1 <- dwstat(errors1, sigma1, modelmatrix1, fittedvalues1, max.lag=8)
	plot(cbind(dw1$r,dw1$p),main= 'Durbin Watson for AR1', xlab= 'Autocorrelation', ylab= 'P value')
	detach("current")
	
# Now try AR2	
	
	ar2 <- gls(bold~hrf+dhrf,correlation=corARMA(p=2))
	auto2 <- coef(ar2$modelStruct$corStruct, F)
	errors2 <- innovations(ar2$resid,auto2)
	current <- as.data.frame(errors2)
	attach(current)
	tsplot(errors2)
	title(paste(dataname,": Innovations(AR2)",sep=""))
	output2 <- acf(errors2, type = "partial", plot = T)
	lm2 <- lm(errors2~1)
	errors2 <- residuals(lm2)
	sigma2 <- summary(lm2)$sigma
	modelmatrix2 <- eval(model.matrix(lm2))
	fittedvalues2 <- fitted.values(lm2)
	dw2 <- dwstat(errors2, sigma2, modelmatrix2, fittedvalues2, max.lag=8)
	plot(cbind(dw2$r,dw2$p),main= 'Durbin Watson for AR2', xlab= 'Autocorrelation', ylab= 'P value')
	detach("current")
	
# Now try AR3	
	
	ar3 <- gls(bold~hrf+dhrf,correlation=corARMA(p=3))
	auto3 <- coef(ar3$modelStruct$corStruct, F)
	errors3 <- innovations(ar3$resid,auto3)
	current <- as.data.frame(errors3)
	attach(current)
	tsplot(errors3)
	title(paste(dataname,": Innovations(AR3)",sep=""))
	output3 <- acf(errors3, type = "partial", plot = T)
	lm3 <- lm(errors3~1)
	errors3 <- residuals(lm3)
	sigma3 <- summary(lm3)$sigma
	modelmatrix3 <- eval(model.matrix(lm3))
	fittedvalues3 <- fitted.values(lm3)
	dw3 <- dwstat(errors3, sigma3, modelmatrix3, fittedvalues3, max.lag=8)
	plot(cbind(dw3$r,dw3$p),main= 'Durbin Watson for AR3', xlab= 'Autocorrelation', ylab= 'P value')
	detach("current")
	
#	print(summary(ols))
#	print(coefficients(ols))
#	print(summary(ar1))
#	print(coefficients(ar1))
#	print(summary(ar2))
#	print(coefficients(ar2))
#	print(summary(ar3))
#	print(coefficients(ar3))

	return(
		AR0=summary(ols), dw0=dw0, Cols=coefficients(ols),
		AR1=summary(ar1), dw1=dw1, Car1=coefficients(ar1),
		AR2=summary(ar2), dw2=dw2, Car2=coefficients(ar2),
		AR3=summary(ar3), dw3=dw3, Car3=coefficients(ar3)
	)
}
}}}