Análise 1st e 2nd level

Teste para um grupo:

Para testar se esse \(\beta\), basta testar com teste-T se o \(\alpha>0\) considerando o modelo de $ Y = + _t$ (\(\epsilon\) com média zero):

\[\begin{align} Y = \begin{bmatrix} \\ \beta_1 \\ \beta_2 \\ ... \\ \beta_n \end{bmatrix} \end{align}\]

e

\[\begin{align} X = \begin{bmatrix} \\ 1 \\ 1 \\ ... \\ 1 \end{bmatrix} \end{align}\]

Isso é equivalente a um teste T para uma amostra.

PRIMEIRO método: Comparando dois grupos:

Por exemplo sejam as matrizes:

\[\begin{align} X = \begin{bmatrix} Controle \\ \beta_1 \\ \beta_2 \\ ... \\ \beta_n \end{bmatrix} \end{align}\]

e pacientes:

\[\begin{align} X = \begin{bmatrix} Pacientes \\ \beta_1 \\ \beta_2 \\ ... \\ \beta_n \end{bmatrix} \end{align}\]

O modelo de GLM utilizado é:

\(Y_i = \beta_1.C_i + \beta_2.P_i + \epsilon_t\) note que não há intercepto \(\alpha\)
e: \(Y_i\) é o coeficiente de ativação do indivíduo:
\(i=\) {1, 2, …, \(N_c\), \(N_c\), …, \((N_c+N_p)\)
\(C_i=1\) se o indivíduo \(i\) é controle e zero caso contrário \(P_i=1\) se o indivíduo é paciente e zero caso contrário

\[ Y= \left(\begin{array}{cc} \beta_1^c \\ \beta_2^c \\ ... \\ \beta_N^p \end{array}\right), X= \left(\begin{array} 1 1 & 0\\ 1 & 0 \\ 1 & 0 \\ ... & ...\\ 0 & 1 \\ 0 & 1 \\ 0 & 1 \\ \end{array}\right) \]

Considerando a base a qual os 37 primeiros indivíduos são controles saudáveis e os demais são pacientes com Parkinson.

Implementando no R:

#Leitura de biblioteca
require(RNifti)

## Carregando pacotes exigidos: RNifti

#Leitura dos dados
betas = readNifti("./data/grupo37C-P.nii")
dim(betas)

## [1] 45 54 45 92

# imagem axial do décimo individuo:
#imagem axial (z=20) do decimo individuo na fatia 20
image(betas[,, 20, 10])

Para um único voxel 20, 20, 20:

Y = betas[20, 20, 20,]
X = matrix(0, 92, 2)
# Como 1: 37 são controle e do 38:92 são pacientes:

# Variável dummy para controles
X[1:37, 1] = 1

# Variável dummy para pacientes:
X[38:92, 2] = 1

Ajustando o modelo - GLM:

#Ajustar o GLM
#o -1 na formula serve para tirar o intercepto
modelo = lm(Y~-1+X[, 1] + X[, 2])
summary(modelo)

## 
## Call:
## lm(formula = Y ~ -1 + X[, 1] + X[, 2])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50669 -0.15868 -0.03905  0.09441  1.15112 
## 
## Coefficients:
##        Estimate Std. Error t value Pr(>|t|)    
## X[, 1]  0.48435    0.04685   10.34   <2e-16 ***
## X[, 2]  0.44014    0.03843   11.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.285 on 90 degrees of freedom
## Multiple R-squared:  0.7256, Adjusted R-squared:  0.7195 
## F-statistic:   119 on 2 and 90 DF,  p-value: < 2.2e-16

Se olhar o \(\beta_1\) são as variáveis do grupo controle e \(\beta_2\) do grupo de pacientes.

Analisando o teste de hipóteses e a ativação dos respectivos controles. Para ativação \(beta_1 > 0\). Olhando o resultado vemos que o p-valor é muito pequeno, portanto, os controles ativam nesse voxel. O mesmo vale para o grupo de pacientes.

Para entender se o controle ativa MAIS que o de pacientes: \(\beta_1 > \beta_2\). No R esse teste não é tão simples de se verificar.

Outro método: Grupos de referência

\(Y_i = \alpha + \beta_2.P_i + \epsilon_t\) \(P_i=1\) se o indivíduo é paciente e zero caso contrário

Para controles: \(Y_i = \alpha + \epsilon_i\)

Para pacientes: \(Y_i= \alpha+ \beta + \epsilon_i\) \(\beta\) é o incremento na média em relação aos controles.

Neste caso:

\[ Y= \left(\begin{array}{cc} \beta_1 \\ \beta_2 \\ ... \\ \beta_{(N_c+N_p)} \end{array}\right), X= \left(\begin{array} 1 1 & 0\\ 1 & 0 \\ 1 & 0 \\ ... & ...\\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ \end{array}\right) \]

Implementando para o mesmo voxel 20,20,20:

### METODO GRUPO DE REFERENCIA
#Neste caso sera o controle

#Observacao: os 37 primeiros indiv
#sao controles saudáveis e os demais
#sao pacientes com Parkinson
#GLM

Y = betas[20, 20, 20,]
X = matrix(0, 92, 2)
#Intercepto
X[, 1] = 1
#Variavel dummy para Pacientes
X[38:92, 2] = 1

#Ajustar o GLM
modelo = lm(Y~X[, 2])
summary(modelo)

## 
## Call:
## lm(formula = Y ~ X[, 2])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.50669 -0.15868 -0.03905  0.09441  1.15112 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.48435    0.04685   10.34   <2e-16 ***
## X[, 2]      -0.04421    0.06060   -0.73    0.468    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.285 on 90 degrees of freedom
## Multiple R-squared:  0.00588,    Adjusted R-squared:  -0.005166 
## F-statistic: 0.5323 on 1 and 90 DF,  p-value: 0.4675

Analisando especificamente para esse voxel, como \(\beta<0\) é como se o grupo de pacientes possui uma ligeira ativação menor que o grupo de controle, Mas ao verificar o p-valor, pode-se verificar que não é significante.

Para análise dos mapas, analisamos mais o t-valor.

Fazendo a varredura para todos os voxeis intracranianos (análise de grupo de referência):

#GLM em todos os voxels intracranianos
MAPAT = array(0, c(45, 54, 45, 1))

for(xi in 1:45){
  for(yi in 1:54){
    for(zi in 1:45){ 
     if(betas[xi, yi, zi, 1] != 0){
      Y = betas[xi, yi, zi, ]
      X = matrix(0, 92, 2)
      #Intercepto
      X[, 1] = 1
      #Variavel dummy para Pacientes
      X[38:92, 2] = 1

      #Ajustar o GLM
      modelo = lm(Y~X[, 2])
      
      #Estatistica T do paciente vs controle
      MAPAT[xi, yi, zi, 1] = summary(modelo)$coef[2, 3]
    }#fecha if
}}}#fecha o for do xi,yi,zi

#Salva os mapas em arquivos no formato Analyze (IMG/HDR)
writeNifti(MAPAT,
           "./data/MapaControlevsPaciente.nii.gz",
           datatype="float")