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)
)
}
