Nuclear Magnetic Resonance Spectroscopy Measurement of Gamma-Amino-Butyric Acid Is Degraded by Head Movement

Motivation

Nuclear Magnetic Resonance Spectroscopy (NMRS) is an important tool for measuring neurochemical levels in human subjects in clinical research. However, its measurement is heavily impacted by subject’s head movement. Here, we want to explore how specifically head movement can impact NMRS measurements.

A brief introduction to data analyzed

Movement data

The size of the raw data is 2.6GB. The head movement data is acquired by video recording of visual markers (Figure A) affixed onto the head of human subjects, with video sampling at 5 Hz (example image shown in Figure B). The position of the markers are extracted by image processing protocol written in MATLAB (https://github.com/sapphire008/MATLAB/tree/master/image_reg3), using techniques including template matching (Modified Hausdorff Displacement) and kmeans clustering. Displacement and Speed (Figure D, showing both vertical and horizontal movements) are obtained from the relative position and the first derivative of position of the centroid of the square images, respectively. This time series is then low-pass filtered with a Butterworth filter at 0.25 hZ cutoff. In the following analyses, however, Displacement and Speed refers to a feature engineered parameter, which is the Euclidean Displacement of the two dimensions over each time point and then taken the root-mean-square of the time series to obtain a single number.

Neurochemical measurements

This is a set of target parameters, including various NMRS measurement of neurochemicals in the brain. The size of the raw data is 10GB.

library(ggplot2)
library(bbmle)
library(plyr)

# Loading existing data
load("D:/Edward/Documents/Assignments/Imaging Research Center/GABA Movement Second Submission/Analysis/data/data_results_11202014.RData")

Building a maximum likelihood model for each neurochemical signals against movement parameters

Qeustion: Do the Displacement and/or speed of head movement negatively impact neurochemical signal measurement?

# Use maximum likelihood to do simple linear regression
# Custom fitting function
LL <- function(params){
  beta = matrix(NA, nrow = length(params) - 2, ncol = 1)
  beta[,1] = params[1:(length(params)-2)]
  mu       = params[[length(params)-1]]
  sigma    = params[[length(params)]]
  minusll  = -sum(suppressWarnings(dnorm(Y - X %*% beta, 0, sigma, log=T))) # Minus of log likelihood
  return(minusll)
}
residuals.mle<- function(params) return(Y -X %*%as.matrix(params[1:(length(params)-2)]))
predict.mle <- function(params) return(X %*%as.matrix(params[1:(length(params)-2)]))

# GABA: against both Displacement and speed of movement
X <- model.matrix(lm(GABA.Cr ~ D + S, data = data.clean.GABA.DS))
Y <- data.clean.GABA.DS$GABA.Cr
parnames(LL) <- c('Intercept', 'D', 'S', 'mu', 'sigma')
mle.GABA.DS <- mle2(LL, start = c(Intercept = 0.1, D = 0.2, S = 0.1, mu = 0, sigma = 1),
     vecpar = TRUE, parnames = c('Intercept', 'D', 'S', 'mu', 'sigma'))
summary(mle.GABA.DS)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, D = 0.2, S = 0.1,
    mu = 0, sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "D", "S", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.6877e-01  6.6436e-03 25.4031 < 2.2e-16 ***
D         -9.5368e-03  4.7700e-03 -1.9993  0.045574 *
S         -1.8287e-01  6.1756e-02 -2.9612  0.003064 **
mu         0.0000e+00  3.0317e-15  0.0000  1.000000
sigma      1.1007e-02  1.4784e-03  7.4454 9.667e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -173.2451 

# GABA: against Displacement of movement only
X <- model.matrix(lm(GABA.Cr ~ D, data = data.clean.GABA.D))
Y <- data.clean.GABA.D$GABA.Cr
parnames(LL) <- c('Intercept', 'D', 'mu', 'sigma')
mle.GABA.D <- mle2(LL, start = c(Intercept = 0.1, D = 0.2, mu = 0, sigma = 1),
                   vecpar = TRUE, parnames=c('Intercept', 'D', 'mu', 'sigma'))
summary(mle.GABA.D)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, D = 0.2, mu = 0,
    sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "D", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.5577e-01  5.5588e-03 28.0222 < 2.2e-16 ***
D         -1.0391e-02  5.1171e-03 -2.0305    0.0423 *
mu         0.0000e+00  9.6376e-16  0.0000    1.0000
sigma      1.2387e-02  1.6345e-03  7.5788 3.487e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -172.5711 

# GABA: against speed of movement only
X <- model.matrix(lm(GABA.Cr ~ S, data = data.clean.GABA.S))
Y <- data.clean.GABA.S$GABA.Cr
parnames(LL) <- c('Intercept', 'S', 'mu', 'sigma')
mle.GABA.S <- mle2(LL, start = c(Intercept = 0.1, S = 0.2, mu = 0, sigma = 1),
                   vecpar = TRUE, parnames=c('Intercept', 'S', 'mu', 'sigma'))
summary(mle.GABA.S)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, S = 0.2, mu = 0,
    sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "S", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.6043e-01  5.1357e-03 31.2385 < 2.2e-16 ***
S         -1.8650e-01  6.0571e-02 -3.0791  0.002076 **
mu         0.0000e+00  2.8010e-14  0.0000  1.000000
sigma      1.1508e-02  1.5228e-03  7.5576 4.105e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -176.9467 

# NAA
X <- model.matrix(lm(NAA.Cr ~ D + S, data = data.clean.NAA.DS))
Y <- data.clean.NAA.DS$NAA.Cr
parnames(LL) <- c('Intercept', 'D', 'S', 'mu', 'sigma')
mle.NAA.DS <- mle2(LL, start = c( Intercept = 0.1, D = 0.2, S = 0.1,  mu = 0, sigma = 1),
                   vecpar = TRUE, parnames=c('Intercept', 'D', 'S', 'mu', 'sigma'))
summary(mle.NAA.DS)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, D = 0.2, S = 0.1,
    mu = 0, sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "D", "S", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value   Pr(z)
Intercept  1.8018e+00  5.9850e-02 30.1058 < 2e-16 ***
D         -5.0059e-02  4.2585e-02 -1.1755  0.2398
S          5.9440e-03  5.7699e-01  0.0103  0.9918
mu         0.0000e+00  1.3518e-11  0.0000  1.0000
sigma      9.8290e-02  1.3136e-02  7.4824 7.3e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -50.45492 

X <- model.matrix(lm(NAA.Cr ~ D, data = data.clean.NAA.D))
Y <- data.clean.NAA.D$NAA.Cr
parnames(LL) <- c('Intercept', 'D', 'mu', 'sigma')
mle.NAA.D <- mle2(LL, start = c(Intercept = 0.1, D = 0.2, mu = 0, sigma = 1),
                   vecpar = TRUE, parnames=c('Intercept', 'D', 'mu', 'sigma'))
summary(mle.NAA.D)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, D = 0.2, mu = 0,
    sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "D", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.7875e+00  4.5878e-02 38.9626 < 2.2e-16 ***
D         -2.8772e-02  4.2114e-02 -0.6832    0.4945
mu         0.0000e+00  6.1171e-14  0.0000    1.0000
sigma      1.0190e-01  1.3382e-02  7.6149 2.639e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -50.16316 

X <- model.matrix(lm(NAA.Cr ~ S, data = data.clean.NAA.S))
Y <- data.clean.NAA.S$NAA.Cr
parnames(LL) <- c('Intercept', 'S', 'mu', 'sigma')
mle.NAA.S <- mle2(LL, start = c(Intercept = 0.1, S = 0.2, mu = 0, sigma = 1),
                  vecpar = TRUE, parnames=c('Intercept', 'S', 'mu', 'sigma'))
summary(mle.NAA.S)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, S = 0.2, mu = 0,
    sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "S", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.7679e+00  4.5673e-02 38.7084 < 2.2e-16 ***
S         -1.2854e-01  5.4529e-01 -0.2357    0.8136
mu         0.0000e+00  7.7801e-13  0.0000    1.0000
sigma      9.9816e-02  1.3108e-02  7.6151 2.636e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -51.36223 

# Glx
X <- model.matrix(lm(Glx.Cr ~ D + S, data = data.clean.Glx.DS))
Y <- data.clean.Glx.DS$Glx.Cr
parnames(LL) <- c('Intercept', 'D', 'S', 'mu', 'sigma')
mle.Glx.DS <- mle2(LL, start = c( Intercept = 0.1, D = 0.2, S = 0.1,  mu = 0, sigma = 1),
                   vecpar = TRUE, parnames=c('Intercept', 'D', 'S', 'mu', 'sigma'))
summary(mle.Glx.DS)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, D = 0.2, S = 0.1,
    mu = 0, sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "D", "S", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.3020e-01  7.3627e-03 17.6838 < 2.2e-16 ***
D         -4.0622e-03  5.2388e-03 -0.7754    0.4381
S          3.6026e-02  7.0981e-02  0.5075    0.6118
mu         0.0000e+00  1.1135e-13  0.0000    1.0000
sigma      1.2092e-02  1.6256e-03  7.4382 1.021e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -168.0179 

X <- model.matrix(lm(Glx.Cr ~ D, data = data.clean.Glx.D))
Y <- data.clean.Glx.D$Glx.Cr
parnames(LL) <- c('Intercept', 'D', 'mu', 'sigma')
mle.Glx.D <- mle2(LL, start = c(Intercept = 0.1, D = 0.2, mu = 0, sigma = 1),
                  vecpar = TRUE, parnames=c('Intercept', 'D', 'mu', 'sigma'))
summary(mle.Glx.D)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, D = 0.2, mu = 0,
    sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "D", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.3287e-01  5.5395e-03 23.9864 < 2.2e-16 ***
D         -4.2456e-03  5.0850e-03 -0.8349    0.4038
mu         0.0000e+00  2.6552e-16  0.0000    1.0000
sigma      1.2304e-02  1.6204e-03  7.5930 3.127e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -172.8912 

X <- model.matrix(lm(Glx.Cr ~ S, data = data.clean.Glx.S))
Y <- data.clean.Glx.S$Glx.Cr
parnames(LL) <- c('Intercept', 'S', 'mu', 'sigma')
mle.Glx.S <- mle2(LL, start = c(Intercept = 0.1, S = 0.2, mu = 0, sigma = 1),
                  vecpar = TRUE, parnames=c('Intercept', 'S', 'mu', 'sigma'))
summary(mle.Glx.S)

Maximum likelihood estimation

Call:
mle2(minuslogl = LL, start = c(Intercept = 0.1, S = 0.2, mu = 0,
    sigma = 1), vecpar = TRUE, parnames = c("Intercept",
    "S", "mu", "sigma"))

Coefficients:
             Estimate  Std. Error z value     Pr(z)
Intercept  1.3132e-01  6.4165e-03 20.4657 < 2.2e-16 ***
S         -4.7113e-02  7.6607e-02 -0.6150    0.5386
mu         0.0000e+00  1.0891e-14  0.0000    1.0000
sigma      1.4023e-02  1.8475e-03  7.5902 3.195e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: -165.3215 

Conclusion: GABA seems to be very sensitive to movement, whereas the other two more abundant chemical does not seem to be affected.

Plot measurements of neurochemical signals in primary visual cortex against movement parameter and fit a line

# GABA
lm.GABA.D <- lm(GABA.Cr ~ D, data.clean.GABA.D)
lm.GABA.D$coefficients <- coef(mle.GABA.D)[1:2]
names(lm.GABA.D$coefficients) <- c("(Intercept)","D")
grid <- with(data.clean.GABA.D, expand.grid(D = seq(min(D), max(D), length=20)))
grid$GABA.Cr <- stats::predict(lm.GABA.D, newdata=grid)
err <- stats::predict(lm.GABA.D, newdata=grid, se = TRUE)
grid$ucl <- err$fit + 1.96 * err$se.fit
grid$lcl <- err$fit - 1.96 * err$se.fit
p = ggplot(data.clean.GABA.D, aes(x=D, y=GABA.Cr)) +
  geom_point(shape=1,size=6) + #geom_smooth(method=lm)+
  geom_smooth(aes(ymin = lcl, ymax = ucl), data = grid, stat="identity")+
  xlab(expression("RMS of Displacement " * Delta ~ " (mm)"))+
  ylab("V1 GABA/Cr")
# L = paste("R^2 ~ \"=\" ~", round(summary(lm(GABA.Cr ~ D, data.clean.GABA.D))$adj.r.squared,digits=4))
L = paste("-2 ~ logL ~ \"=\" ~", round(as.numeric(logLik(mle.GABA.D))*-2,digits=2))
p + geom_text(aes(x=1.5,y = 0.17, label = L),size=6,parse=T)+
  theme(axis.text.x = element_text(size=18),axis.title.x = element_text(size=20),
        axis.text.y = element_text(size=18),axis.title.y = element_text(size=20))

lm.GABA.S <- lm(GABA.Cr ~ S, data.clean.GABA.S)
lm.GABA.S$coefficients <- coef(mle.GABA.S)[1:2]
names(lm.GABA.S$coefficients) <- c("(Intercept)","S")
grid <- with(data.clean.GABA.S, expand.grid(S = seq(min(S), max(S), length=20)))
grid$GABA.Cr <- stats::predict(lm.GABA.S, newdata=grid)
err <- stats::predict(lm.GABA.S, newdata=grid, se = TRUE)
grid$ucl <- err$fit + 1.96 * err$se.fit
grid$lcl <- err$fit - 1.96 * err$se.fit
p = ggplot(data.clean.GABA.S, aes(x=S, y=GABA.Cr)) +
  geom_point(shape=1,size=6) +#geom_smooth(method=lm)+
  geom_smooth(aes(ymin = lcl, ymax = ucl), data = grid, stat="identity")+
  xlab(expression("RMS of Speed " * Sigma ~ " (mm/s)"))+
  ylab("V1 GABA/Cr")
L = paste("-2 ~ logL ~ \"=\" ~", round(as.numeric(logLik(mle.GABA.S))*-2,digits=2))
p + geom_text(aes(x=0.14,y = 0.17, label = L),size=6,parse=T)+
  theme(axis.text.x = element_text(size=18),axis.title.x = element_text(size=20),
        axis.text.y = element_text(size=18),axis.title.y = element_text(size=20))

# NAA
lm.NAA.D <- lm(NAA.Cr ~ D, data.clean.NAA.D)
lm.NAA.D$coefficients <- coef(mle.NAA.D)[1:2]
names(lm.NAA.D$coefficients) <- c("(Intercept)","D")
grid <- with(data.clean.NAA.D, expand.grid(D = seq(min(D), max(D), length=20)))
grid$NAA.Cr <- stats::predict(lm.NAA.D, newdata=grid)
err <- stats::predict(lm.NAA.D, newdata=grid, se = TRUE)
grid$ucl <- err$fit + 1.96 * err$se.fit
grid$lcl <- err$fit - 1.96 * err$se.fit
p = ggplot(data.clean.NAA.D, aes(x=D, y=NAA.Cr)) +
  geom_point(shape=1,size=6) + #geom_smooth(method=lm)+
  geom_smooth(aes(ymin = lcl, ymax = ucl), data = grid, stat="identity")+
  xlab(expression("RMS of Displacement " * Delta ~ " (mm)"))+
  ylab("V1 NAA/Cr")
# L = paste("R^2 ~ \"=\" ~", round(summary(lm(NAA.Cr ~ D, data.clean.NAA.D))$adj.r.squared,digits=4))
L = paste("-2 ~ logL ~ \"=\" ~", round(as.numeric(logLik(mle.NAA.D))*-2,digits=2))
p + geom_text(aes(x=1.5,y = 2.0, label = L),size=6,parse=T)+
  theme(axis.text.x = element_text(size=18),axis.title.x = element_text(size=20),
        axis.text.y = element_text(size=18),axis.title.y = element_text(size=20))

lm.NAA.S <- lm(NAA.Cr ~ S, data.clean.NAA.S)
lm.NAA.S$coefficients <- coef(mle.NAA.S)[1:2]
names(lm.NAA.S$coefficients) <- c("(Intercept)","S")
grid <- with(data.clean.NAA.S, expand.grid(S = seq(min(S), max(S), length=20)))
grid$NAA.Cr <- stats::predict(lm.NAA.S, newdata=grid)
err <- stats::predict(lm.NAA.S, newdata=grid, se = TRUE)
grid$ucl <- err$fit + 1.96 * err$se.fit
grid$lcl <- err$fit - 1.96 * err$se.fit
p = ggplot(data.clean.NAA.S, aes(x=S, y=NAA.Cr)) +
  geom_point(shape=1,size=6) +#geom_smooth(method=lm)+
  geom_smooth(aes(ymin = lcl, ymax = ucl), data = grid, stat="identity")+
  xlab(expression("RMS of Speed " * Sigma ~ " (mm/s)"))+
  ylab("V1 NAA/Cr")
L = paste("-2 ~ logL ~ \"=\" ~", round(as.numeric(logLik(mle.NAA.S))*-2,digits=2))
p + geom_text(aes(x=0.13,y = 2.0, label = L),size=6,parse=T)+
  theme(axis.text.x = element_text(size=18),axis.title.x = element_text(size=20),
        axis.text.y = element_text(size=18),axis.title.y = element_text(size=20))

# Glx
lm.Glx.D <- lm(Glx.Cr ~ D, data.clean.Glx.D)
lm.Glx.D$coefficients <- coef(mle.Glx.D)[1:2]
names(lm.Glx.D$coefficients) <- c("(Intercept)","D")
grid <- with(data.clean.Glx.D, expand.grid(D = seq(min(D), max(D), length=20)))
grid$Glx.Cr <- stats::predict(lm.Glx.D, newdata=grid)
err <- stats::predict(lm.Glx.D, newdata=grid, se = TRUE)
grid$ucl <- err$fit + 1.96 * err$se.fit
grid$lcl <- err$fit - 1.96 * err$se.fit
p = ggplot(data.clean.Glx.D, aes(x=D, y=Glx.Cr)) +
  geom_point(shape=1,size=6) + #geom_smooth(method=lm)+
  geom_smooth(aes(ymin = lcl, ymax = ucl), data = grid, stat="identity")+
  xlab(expression("RMS of Displacement " * Delta ~ " (mm)"))+
  ylab("V1 Glx/Cr")
# L = paste("R^2 ~ \"=\" ~", round(summary(lm(Glx.Cr ~ D, data.clean.Glx.D))$adj.r.squared,digits=4))
L = paste("-2 ~ logL ~ \"=\" ~", round(as.numeric(logLik(mle.Glx.D))*-2,digits=2))
p + geom_text(aes(x=1.5,y = 0.16, label = L),size=6,parse=T)+
  theme(axis.text.x = element_text(size=18),axis.title.x = element_text(size=20),
        axis.text.y = element_text(size=18),axis.title.y = element_text(size=20))

lm.Glx.S <- lm(Glx.Cr ~ S, data.clean.Glx.S)
lm.Glx.S$coefficients <- coef(mle.Glx.S)[1:2]
names(lm.Glx.S$coefficients) <- c("(Intercept)","S")
grid <- with(data.clean.Glx.S, expand.grid(S = seq(min(S), max(S), length=20)))
grid$Glx.Cr <- stats::predict(lm.Glx.S, newdata=grid)
err <- stats::predict(lm.Glx.S, newdata=grid, se = TRUE)
grid$ucl <- err$fit + 1.96 * err$se.fit
grid$lcl <- err$fit - 1.96 * err$se.fit
p = ggplot(data.clean.Glx.S, aes(x=S, y=Glx.Cr)) +
  geom_point(shape=1,size=6) +#geom_smooth(method=lm)+
  geom_smooth(aes(ymin = lcl, ymax = ucl), data = grid, stat="identity")+
  xlab(expression("RMS of Speed " * Sigma ~ " (mm/s)"))+
  ylab("V1 Glx/Cr")
L = paste("-2 ~ logL ~ \"=\" ~", round(as.numeric(logLik(mle.Glx.S))*-2,digits=2))
p + geom_text(aes(x=0.13,y = 0.16, label = L),size=6,parse=T)+
  theme(axis.text.x = element_text(size=18),axis.title.x = element_text(size=20),
        axis.text.y = element_text(size=18),axis.title.y = element_text(size=20))

Is it true that the longer the subject stay in the scanner, the worse the signal?

library(nlme)
load("D:/Edward/Documents/Assignments/Imaging Research Center/GABA Movement Fourth Submission/analysis/data/data_results_10302016.RData")

lme.GABA.Cr = lme(GABA.Cr ~ Region.Run.Order, data=data.homogeneous, random=~1|Subject)
summary(lme.GABA.Cr)

Linear mixed-effects model fit by REML
 Data: data.homogeneous 

Random effects:
 Formula: ~1 | Subject
        (Intercept)    Residual
StdDev: 0.009674584 0.005029841

Fixed effects: GABA.Cr ~ Region.Run.Order 

 Correlation:
                  (Intr) R.R.O2
Region.Run.Order2 -0.326
Region.Run.Order3 -0.326  0.500

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-1.6856021 -0.7082410  0.1739974  0.5781646  1.6793817

Number of Observations: 45
Number of Groups: 15 

lme.Glx.Cr = lme(Glx.Cr ~ Region.Run.Order, data=data.homogeneous, random=~1|Subject)
summary(lme.Glx.Cr)

Linear mixed-effects model fit by REML
 Data: data.homogeneous 

Random effects:
 Formula: ~1 | Subject
        (Intercept)    Residual
StdDev: 0.006411252 0.006937828

Fixed effects: Glx.Cr ~ Region.Run.Order 

 Correlation:
                  (Intr) R.R.O2
Region.Run.Order2 -0.519
Region.Run.Order3 -0.519  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-1.80202593 -0.43097103  0.04998176  0.60270971  2.22699679

Number of Observations: 45
Number of Groups: 15 

lme.NAA.Cr = lme(NAA.Cr ~ Region.Run.Order, data=data.homogeneous, random=~1|Subject)
summary(lme.NAA.Cr)

Linear mixed-effects model fit by REML
 Data: data.homogeneous 

Random effects:
 Formula: ~1 | Subject
        (Intercept)   Residual
StdDev:   0.1156031 0.02517245

Fixed effects: NAA.Cr ~ Region.Run.Order 

 Correlation:
                  (Intr) R.R.O2
Region.Run.Order2 -0.15
Region.Run.Order3 -0.15   0.50

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-1.73280721 -0.45729043  0.08660261  0.59497552  1.88347078

Number of Observations: 45
Number of Groups: 15 

Conclusion: Yes. The longer the subjected stayed in the scanner (i.e. the longer time it takes to complete the measurement), the worse the GABA and Glx signal, but not NAA signal.