SCodebits - Methods
location: SCodebits

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

None: SCodebits (last edited 2013-03-08 10:28:25 by localhost)