User Tools

Site Tools


neuroimagen:altdti

experimento DTI

pruebas

Siguiendo los pasos de la pipeline 3.0 procesamos las imagenes DTI en intentamos hallar una red que parta de una parte especifica del cortex.

El archivo dti_track.seed se fabrica con 42 regiones distintas,

grep ctx /usr/local/freesurfer/FreeSurferColorLUT.txt | grep 'parietal\|frontal' | awk {'print $1'} > /nas/data/facehbi/dti_track.seed 

Estas regiones son,

1003 ctx-lh-caudalmiddlefrontal
1008 ctx-lh-inferiorparietal
1012 ctx-lh-lateralorbitofrontal
1014 ctx-lh-medialorbitofrontal
1027 ctx-lh-rostralmiddlefrontal
1028 ctx-lh-superiorfrontal
1029 ctx-lh-superiorparietal
1032 ctx-lh-frontalpole
2003 ctx-rh-caudalmiddlefrontal
2008 ctx-rh-inferiorparietal
2012 ctx-rh-lateralorbitofrontal
2014 ctx-rh-medialorbitofrontal
2027 ctx-rh-rostralmiddlefrontal
2028 ctx-rh-superiorfrontal
2029 ctx-rh-superiorparietal
2032 ctx-rh-frontalpole
1106 ctx-lh-G_frontal_inf-Opercular_part
1107 ctx-lh-G_frontal_inf-Orbital_part
1108 ctx-lh-G_frontal_inf-Triangular_part
1109 ctx-lh-G_frontal_middle
1110 ctx-lh-G_frontal_superior
1122 ctx-lh-G_parietal_inferior-Angular_part
1123 ctx-lh-G_parietal_inferior-Supramarginal_part
1124 ctx-lh-G_parietal_superior
1154 ctx-lh-S_frontal_inferior
1155 ctx-lh-S_frontal_middle
1156 ctx-lh-S_frontal_superior
1159 ctx-lh-S_intraparietal-and_Parietal_transverse
1177 ctx-lh-S_subparietal
2106 ctx-rh-G_frontal_inf-Opercular_part
2107 ctx-rh-G_frontal_inf-Orbital_part
2108 ctx-rh-G_frontal_inf-Triangular_part
2109 ctx-rh-G_frontal_middle
2110 ctx-rh-G_frontal_superior
2122 ctx-rh-G_parietal_inferior-Angular_part
2123 ctx-rh-G_parietal_inferior-Supramarginal_part
2124 ctx-rh-G_parietal_superior
2154 ctx-rh-S_frontal_inferior
2155 ctx-rh-S_frontal_middle
2156 ctx-rh-S_frontal_superior
2159 ctx-rh-S_intraparietal-and_Parietal_transverse
2177 ctx-rh-S_subparietal

El resultado de probtrackx es una red bastante extensa,

Ahora voy a quitar los giros y demas y quedarme solo con las 16 primeras lineas.

$ cat dti_track.seed 
1003
1008
1012
1014
1027
1028
1029
1032
2003
2008
2012
2014
2027
2028
2029
2032

La red que obtengo es ahora bastante similar,

Aplicando esta red como mascara (25%) saco los datos de FA en la red y los puedo comparar con otros datos . Ejemplo, SUVR en la misma visita,

Aislando el cortex

Voy a intentar aislar todas las regiones del cortex para hacer una lista de lo que puedo parcelar

$ grep ctx /usr/local/freesurfer/FreeSurferColorLUT.txt | grep -v "G\|S\|_\|nknown\|#" | awk {'if($1<3000) print $1,$2'} | grep lh | sed 's/-lh//' | awk {'print $2'}

y me quedo con los nombres utiles para que un neurologo me edite la lista,

ctx-bankssts
ctx-caudalanteriorcingulate
ctx-caudalmiddlefrontal
ctx-corpuscallosum
ctx-cuneus
ctx-entorhinal
ctx-fusiform
ctx-inferiorparietal
ctx-inferiortemporal
ctx-isthmuscingulate
ctx-lateraloccipital
ctx-lateralorbitofrontal
ctx-lingual
ctx-medialorbitofrontal
ctx-middletemporal
ctx-parahippocampal
ctx-paracentral
ctx-parsopercularis
ctx-parsorbitalis
ctx-parstriangularis
ctx-pericalcarine
ctx-postcentral
ctx-posteriorcingulate
ctx-precentral
ctx-precuneus
ctx-rostralanteriorcingulate
ctx-rostralmiddlefrontal
ctx-superiorfrontal
ctx-superiorparietal
ctx-superiortemporal
ctx-supramarginal
ctx-frontalpole
ctx-temporalpole
ctx-transversetemporal
ctx-insula

Mapa FP MB

Las regiones escogidas son,

ctx-caudalmiddlefrontal
ctx-inferiorparietal
ctx-middletemporal
ctx-parsopercularis
ctx-parstriangularis
ctx-postcentral
ctx-precentral
ctx-superiorfrontal
ctx-superiorparietal
ctx-superiortemporal
ctx-supramarginal

Voy a quedarme con lo que necesito ahora,

$ grep ctx /usr/local/freesurfer/FreeSurferColorLUT.txt | grep -v "G\|S\|_\|nknown\|#" | awk {'if($1<3000) print $1,$2'} > tofp.txt
$ cat tocut.txt
caudalmiddlefrontal
inferiorparietal
middletemporal
parsopercularis
parstriangularis
postcentral
precentral
superiorfrontal
superiorparietal
superiortemporal
supramarginal
$ grep -f tocut.txt tofp.txt | awk {'print $1'} > dti_track.seed 
$ cat dti_track.seed
1003
1008
1015
1018
1020
1022
1024
1028
1029
1030
1031
2003
2008
2015
2018
2020
2022
2024
2028
2029
2030
2031

con estas seeds ya puedo correr el experimento.

$ dti_track.pl facehbi
... 3dias ...
$ cd /nas/data/facehbi
$ for x in `ls -d working/*_probtrack_out`; do mv $x `echo $x | sed 's/out/FPCustom/'`;done
$ dti_metrics_alt.pl -path FPCustom facehbi
 
$ dti_track.pl v2MriPet
... 3dias ...
$ cd /nas/data/v2MriPet
$ for x in `ls -d working/*_probtrack_out`; do mv $x `echo $x | sed 's/out/FPCustom/'`;done
$ dti_metrics_alt.pl -path FPCustom v2MriPet

Comparando con DMN

Sorprendentemente, aunque estan relacionados, hay bastante diferencia entre la FA medida en la DMN y la medida en la nueva red,

> summary(v1m)
 
Call:
lm(formula = DMN_FA_v1 ~ FPCustom_FA_v1, data = dti_c)
 
Residuals:
      Min        1Q    Median        3Q       Max 
-0.031888 -0.006546  0.001095  0.006806  0.026302 
 
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.06472    0.01262   5.129 8.74e-07 ***
FPCustom_FA_v1  0.89746    0.04285  20.943  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.01099 on 152 degrees of freedom
Multiple R-squared:  0.7426,	Adjusted R-squared:  0.741 
F-statistic: 438.6 on 1 and 152 DF,  p-value: < 2.2e-16

> summary(v2m)
 
Call:
lm(formula = DMN_FA_v2 ~ FPCustom_FA_v2, data = dti_c)
 
Residuals:
      Min        1Q    Median        3Q       Max 
-0.036853 -0.005949  0.000129  0.006172  0.024941 
 
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     0.07965    0.01236   6.443 1.46e-09 ***
FPCustom_FA_v2  0.83691    0.04172  20.058  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.01042 on 152 degrees of freedom
Multiple R-squared:  0.7258,	Adjusted R-squared:  0.724 
F-statistic: 402.3 on 1 and 152 DF,  p-value: < 2.2e-16

Modelo mixto

Nada raro, la relacion con el SUVR va mas o menos igual.

> model.a <- lm(Custom_FA ~ SUVR , data=idata)
> summary(model.a)
 
Call:
lm(formula = Custom_FA ~ SUVR, data = idata)
 
Residuals:
      Min        1Q    Median        3Q       Max 
-0.055969 -0.015556 -0.001559  0.015067  0.055285 
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.308468   0.007587  40.659   <2e-16 ***
SUVR        -0.011063   0.006002  -1.843   0.0662 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.02038 on 306 degrees of freedom
Multiple R-squared:  0.01098,	Adjusted R-squared:  0.00775 
F-statistic: 3.398 on 1 and 306 DF,  p-value: 0.06625
 
> model.c <- lme(Custom_FA ~ SUVR, random = ~ 1| Subject, data=idata)
> summary(model.c)
Linear mixed-effects model fit by REML
 Data: idata 
       AIC       BIC logLik
  -1566.58 -1551.686 787.29
 
Random effects:
 Formula: ~1 | Subject
        (Intercept)  Residual
StdDev:  0.01580869 0.0129124
 
Fixed effects: Custom_FA ~ SUVR 
                  Value   Std.Error  DF  t-value p-value
(Intercept)  0.30835579 0.008911474 153 34.60211  0.0000
SUVR        -0.01097326 0.007035937 153 -1.55960  0.1209
 Correlation: 
     (Intr)
SUVR -0.986
 
Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.10991932 -0.50342918 -0.01151992  0.50847071  1.97546782 
 
Number of Observations: 308
Number of Groups: 154
 
> anova(model.c, model.a)
        Model df       AIC       BIC   logLik   Test  L.Ratio p-value
model.c     1  4 -1566.580 -1551.686 787.2900                        
model.a     2  3 -1500.089 -1488.918 753.0443 1 vs 2 68.49128  <.0001
 
> model.b <- lmer(Custom_FA ~ SUVR + (1| Subject), data=idata)
> summary(model.b)
Linear mixed model fit by REML ['lmerMod']
Formula: Custom_FA ~ SUVR + (1 | Subject)
   Data: idata
 
REML criterion at convergence: -1574.6
 
Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.10992 -0.50343 -0.01152  0.50847  1.97547 
 
Random effects:
 Groups   Name        Variance  Std.Dev.
 Subject  (Intercept) 0.0002499 0.01581 
 Residual             0.0001667 0.01291 
Number of obs: 308, groups:  Subject, 154
 
Fixed effects:
             Estimate Std. Error t value
(Intercept)  0.308356   0.008911   34.60
SUVR        -0.010973   0.007036   -1.56
 
Correlation of Fixed Effects:
     (Intr)
SUVR -0.986
 
> anova(model.b, model.a)
refitting model(s) with ML (instead of REML)
Data: idata
Models:
model.a: Custom_FA ~ SUVR
model.b: Custom_FA ~ SUVR + (1 | Subject)
        Df     AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
model.a  3 -1520.2 -1509 763.08  -1526.2                             
model.b  4 -1585.9 -1571 796.94  -1593.9 67.704      1  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Modelos

Estratificando por APOE

Tenemos los genotipos de APOE en un CSV,

> head updateAPOE_FACEHBI_051218.csv
code_facehbi;Interno;APOE
F079;20120457;e3e3
F103;20121210;e3e3
F080;20121207;e2e2
F097;20130137;e3e3
F018;20130447;e3e4
F096;20130432;e3e4
F002;20131084;e3e3
F113;20131063;e3e4
F027;20131385;e3e4

Limpiamos un poco esto,

> awk -F";" {'print $1,$3'} updateAPOE_FACEHBI_051218.csv | sed 's/F//; s/ /,/; s/code_facehbi/Subject/' | sort -n > facehbi_apoe.csv
> head facehbi_apoe.csv 
Subject,APOE
001,e3e4
002,e3e3
003,e3e3
004,e3e3
005,e2e3
006,e3e4
007,e3e3
008,e3e3
009,e3e3

Lo separamos para estratificar,

> sed 's/e2e2\|e2e3/0/; s/e2e4\|e3e3/1/; s/e3e4\|e4e4/2/' facehbi_apoe.csv > facehbi_apoe_strats.csv
> head facehbi_apoe_strats.csv
Subject,APOE
001,2
002,1
003,1
004,1
005,0
006,2
007,1
008,1
009,1

Esto lo importo en R (mas o menos),

> idatawnp <- read.csv("facehbi_dmn_fa_suvr_np_tmp_v1.csv", sep=";", header=TRUE)
> iapoe <-read.csv("facehbi_apoe_strats.csv", sep=";", header= TRUE)
> idatacustom <- read.csv("facehbi_suvr_dmnfa_customfa_v12.csv", sep=";", header=TRUE)
> idatatmp <-merge(idatawnp, idatacustom, by="Subject")
> idata <-merge(idatatmp, iapoe, by="Subject")
> idata_apoe0 <- idata[idata$APOE == "0",]

Un monton de data aqui pero nos centramos en la visita basal,

> m0 <- lm(idata_apoe0$FPCustom_FA_v1 ~ idata_apoe0$Global_v1.x)
> summary(m0)
 
Call:
lm(formula = idata_apoe0$FPCustom_FA_v1 ~ idata_apoe0$Global_v1.x)
 
Residuals:
      Min        1Q    Median        3Q       Max 
-0.051447 -0.010969  0.001284  0.009302  0.040698 
 
Coefficients:
                        Estimate Std. Error t value Pr(>|t|)  
(Intercept)              0.18674    0.07971   2.343   0.0278 *
idata_apoe0$Global_v1.x  0.09175    0.06785   1.352   0.1889  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.02192 on 24 degrees of freedom
Multiple R-squared:  0.07078,	Adjusted R-squared:  0.03207 
F-statistic: 1.828 on 1 and 24 DF,  p-value: 0.1889

OK. ahi no hay nada pero esto solo era para ver si funcionaba la cosa. Vamos a organizar un poco la tabla,

> colnames(idata)
 [1] "Subject"                        "DMN_FA_v1.x"                    "Global_v1.x"                   
 [4] "Escolaridad"                    "male"                           "LL_Namingtotal_NP"             
 [7] "LL_total_NP"                    "total_NP"                       "Piramides_Y_Palmeras_Palab_FAC"
[10] "Piramides_Y_Plameras_Imag_FAC"  "Kissing_Dancing_Imagenes_FAC"   "Kissing_Dancing_Palabras_FAC"  
[13] "FNAME_TOTAL_FAC"                "Action_Namimg_Libre_FAC"        "Boston_Libre_FAC"              
[16] "Boston_Total_FAC"               "Global_v1.y"                    "Global_v2"                     
[19] "DMN_FA_v1.y"                    "DMN_FA_v2"                      "FPCustom_FA_v1"                
[22] "FPCustom_FA_v2"                 "APOE" 

Con esos nombres me voy a equivocar seguro, asi que

> colnames(idata)[colnames(idata) == "Piramides_Y_Plameras_Imag_FAC"] <- "PyP_i"
> colnames(idata)[colnames(idata) == "Piramides_Y_Palmeras_Palab_FAC"] <- "PyP_p"
> colnames(idata)[colnames(idata) == "Kissing_Dancing_Imagenes_FAC"] <- "KD_i"
> colnames(idata)[colnames(idata) == "Kissing_Dancing_Palabras_FAC"] <- "KD_p"
> colnames(idata)[colnames(idata) == "LL_Namingtotal_NP"] <- "LL_Naming"
> colnames(idata)[colnames(idata) == "LL_total_NP"] <- "LL"
> colnames(idata)[colnames(idata) == "total_NP"] <- "NP"
> colnames(idata)[colnames(idata) == "FNAME_TOTAL_FAC"] <- "FName"
> colnames(idata)[colnames(idata) == "Action_Namimg_Libre_FAC"] <- "ANaming"
> colnames(idata)[colnames(idata) == "Boston_Libre_FAC"] <- "Boston_Libre"
> colnames(idata)[colnames(idata) == "Boston_Total_FAC"] <- "Boston"
> colnames(idata)
 [1] "Subject"        "DMN_FA_v1.x"    "Global_v1.x"    "Escolaridad"    "male"           "LL_Naming"     
 [7] "LL"             "NP"             "PyP_p"          "PyP_i"          "KD_i"           "KD_p"          
[13] "FName"          "ANaming"        "Boston_Libre"   "Boston"         "Global_v1.y"    "Global_v2"     
[19] "DMN_FA_v1.y"    "DMN_FA_v2"      "FPCustom_FA_v1" "FPCustom_FA_v2" "APOE" 

y tambien,

> colnames(idata)[colnames(idata) == "DMN_FA_v1.x"] <- "DMN"
> colnames(idata)[colnames(idata) == "Global_v1.x"] <- "SUVR"
> colnames(idata)[colnames(idata) == "FPCustom_FA_v1"] <- "FPCustom"
> drops <- c("Global_v1.y", "Global_v2", "DMN_FA_v1.y", "DMN_FA_v2", "FPCustom_FA_v2")
> idata[ , !(names(idata) %in% drops)]
> okdata <- idata[ , !(names(idata) %in% drops)]

con lo cual queda mucho mas potable esto,

> colnames(okdata)
 [1] "Subject"      "DMN"          "SUVR"         "Escolaridad"  "male"         "LL_Naming"    "LL"          
 [8] "NP"           "PyP_p"        "PyP_i"        "KD_i"         "KD_p"         "FName"        "ANaming"     
[15] "Boston_Libre" "Boston"       "FPCustom"     "APOE"  

Ahora, el objetivo es estudiar las variables neurosicologicas como funcion de la FA en las redes DTI (DMN y FPCustom), estratificando por APOE y usando como covariables SUVR, Genero, Escolaridad y Edad. m( joer q me falta la edad.

> awk -F";" '{print $2,$8}' faceHBI_matriuREF_14-1-19-1.csv | sed 's/ /;/; s/edat/Edad/; s/subject/Subject/' > edad.csv
> scp -P 20022 edad.csv detritus.fundacioace.org:facehbi/dti_model/
edad.csv                                                      100% 1312   200.6KB/s   00:00  . 

Vale,la meto

> edad_dlc <-read.csv("edad.csv", sep=";", header=TRUE)
> okdata <- merge(okdata, edad_dlc, by = "Subject")
> colnames(okdata)
 [1] "Subject"      "DMN"          "SUVR"         "Escolaridad"  "male"         "LL_Naming"    "LL"          
 [8] "NP"           "PyP_p"        "PyP_i"        "KD_i"         "KD_p"         "FName"        "ANaming"     
[15] "Boston_Libre" "Boston"       "FPCustom"     "APOE"         "Edad" 

Ahora si. Pero sigue siendo una puñeta. Mejor me hago un script que corra a traves de los modelos ;-). Venga, primero a lo bruto, just in case,

library(QuantPsyc)
x<-read.csv("facehbi_dti_np.csv")
Color=c("red","blue")
scan("npvars.names", what = character())->np
scan("nivars.names", what = character())->ni
sink(file = "facehbi_dti_np_models.txt", append = TRUE, type = "output", split = TRUE)
 
for(i in 1:length(np)){
        for(j in 1:length(ni)){
                y.data <- x[c(ni[j], np[i], "male", "Edad", "Escolaridad", "SUVR")]
                y.data <- y.data[complete.cases(y.data),]
                a <- lm( paste ('y.data$', np[i], ' ~ y.data$', ni[j], ' + y.data$male + y.data$Edad + y.data$Escolaridad + y.data$SUVR'))
                writeLines(paste("NP: ", np[i], " NI: ", ni[j]))
                writeLines(paste("R2: ", summary(a)$adj.r.squared, " p-value: ", summary(a)$coef[2,4]))
                beta <- lm.beta(a)
                for(k in 1:length(beta)){
                        writeLines(paste(names(beta[k]), ": ", beta[k]))
                }
                writeLines(paste("-------"))
        }
}
sink()

a ver,

> write.csv(okdata, "facehbi_dti_np.csv")
> source("get_lms.r")
Read 11 items
Read 2 items

Un vistazo a la salida, y effectivamente, no hay nada que hacer aqui. el mas prometedor es este,

NP:  NP  NI:  FPCustom
R2:  0.245411609394026  p-value:  0.0332379182003353

Puaf. Vamos a estratificar aver,

 okdata0 <- okdata[okdata$APOE == "0",]
> write.csv(okdata0, "facehbi_dti_np.csv")
> source("get_lms.r")
Read 11 items
Read 2 items

Poca cosa aqui,

NP:  NP  NI:  DMN
R2:  0.283699960037635  p-value:  0.624010265461206
NP:  FName  NI:  DMN
R2:  0.260050114900356  p-value:  0.141503076825646

Seguimos,

 okdata1 <- okdata[okdata$APOE == "1",]
> write.csv(okdata1, "facehbi_dti_np.csv")
> source("get_lms.r")
Read 11 items
Read 2 items

Nop.

 okdata2 <- okdata[okdata$APOE == "2",]
> write.csv(okdata2, "facehbi_dti_np.csv")
> source("get_lms.r")
Read 11 items
Read 2 items

grrrrrrr…

NP:  LL_Naming  NI:  DMN
R2:  0.338308455461783  p-value:  0.00172399799381871
NP:  NP  NI:  FPCustom
R2:  0.423468586516873  p-value:  0.344571425163692
NP:  KD_p  NI:  FPCustom
R2:  0.300292939478652  p-value:  0.382224137819204
NP:  KD_i  NI:  FPCustom
R2:  0.3278128869582  p-value:  0.511816038951647
NP:  FName  NI:  FPCustom
R2:  0.377022589868983  p-value:  0.457838263496454
NP:  Boston  NI:  FPCustom
R2:  0.305029202245965  p-value:  0.251505671812935
NP:  Boston_Libre  NI:  FPCustom
R2:  0.349670089711  p-value:  0.127585742382962

Odio estas cosas, todo porqueria que hay que mirar.

> m0 <- lm(okdata2$LL_Naming ~ okdata2$DMN + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR +okdata2$male)
> summary(m0)
 
Call:
lm(formula = okdata2$LL_Naming ~ okdata2$DMN + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-1.0747 -0.2115  0.1003  0.2248  0.4808 
 
Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          20.55919    1.58086  13.005 2.56e-14 ***
okdata2$DMN         -13.72502    4.01223  -3.421  0.00172 ** 
okdata2$Escolaridad   0.04295    0.01474   2.913  0.00648 ** 
okdata2$Edad         -0.03229    0.00971  -3.326  0.00222 ** 
okdata2$SUVR          0.16477    0.26430   0.623  0.53743    
okdata2$male          0.02292    0.13024   0.176  0.86142    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.3551 on 32 degrees of freedom
Multiple R-squared:  0.4277,	Adjusted R-squared:  0.3383 
F-statistic: 4.783 on 5 and 32 DF,  p-value: 0.002238
 
> m0 <- lm(okdata2$LL_Naming ~ okdata2$DMN + okdata2$Escolaridad + okdata2$Edad)
> summary(m0)
 
Call:
lm(formula = okdata2$LL_Naming ~ okdata2$DMN + okdata2$Escolaridad + 
    okdata2$Edad)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-1.13733 -0.21729  0.09543  0.23242  0.42227 
 
Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)          20.529979   1.521191  13.496 3.25e-15 ***
okdata2$DMN         -13.286680   3.840520  -3.460  0.00148 ** 
okdata2$Escolaridad   0.041280   0.014099   2.928  0.00605 ** 
okdata2$Edad         -0.030158   0.008949  -3.370  0.00188 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.3469 on 34 degrees of freedom
Multiple R-squared:  0.4198,	Adjusted R-squared:  0.3686 
F-statistic: 8.201 on 3 and 34 DF,  p-value: 0.0003054

Baaaah, pero he escrito papers con menos. Esta es la variable que se llama en la tabla LLNaming_total_NP, a saber que sera.

Seguimos, la segunda asociacion no dice nada,

> m0 <- lm(okdata2$NP ~ okdata2$FPCustom + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR + okdata2$male)
> summary(m0)

Call:
lm(formula = okdata2$NP ~ okdata2$FPCustom + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)

Residuals:
    Min      1Q  Median      3Q     Max 
-26.514  -7.122   1.120   5.021  32.711 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         176.2271    44.0144   4.004 0.000346 ***
okdata2$FPCustom    121.2372   126.3732   0.959 0.344571    
okdata2$Escolaridad   1.2925     0.5144   2.513 0.017215 *  
okdata2$Edad         -0.7723     0.3146  -2.455 0.019701 *  
okdata2$SUVR         -5.1581     9.0583  -0.569 0.573036    
okdata2$male         -7.5469     4.4865  -1.682 0.102280    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.25 on 32 degrees of freedom
Multiple R-squared:  0.5014,	Adjusted R-squared:  0.4235 
F-statistic: 6.435 on 5 and 32 DF,  p-value: 0.0003044

> m0 <- lm(okdata2$NP ~ okdata2$Escolaridad + okdata2$Edad)
> summary(m0)

Call:
lm(formula = okdata2$NP ~ okdata2$Escolaridad + okdata2$Edad)

Residuals:
    Min      1Q  Median      3Q     Max 
-27.939  -8.188  -0.075   6.659  37.356 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         209.0838    23.1923   9.015 1.19e-10 ***
okdata2$Escolaridad   1.4459     0.4692   3.082  0.00400 ** 
okdata2$Edad         -0.9078     0.3010  -3.016  0.00475 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.41 on 35 degrees of freedom
Multiple R-squared:  0.4404,	Adjusted R-squared:  0.4084 
F-statistic: 13.77 on 2 and 35 DF,  p-value: 3.87e-05

la variable NP esta asociada casi exclusivamente, en este grupo a la edad y la escolaridad del sujeto. Lomismo vale para el resto,

> m0 <- lm(okdata2$KD_p ~ okdata2$FPCustom + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR + okdata2$male)
> summary(m0)

Call:
lm(formula = okdata2$KD_p ~ okdata2$FPCustom + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0610 -0.7722  0.1195  0.6817  2.1740 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          59.01916    8.95276   6.592 6.19e-06 ***
okdata2$FPCustom    -23.02432   25.62433  -0.899   0.3822    
okdata2$Escolaridad   0.27985    0.09696   2.886   0.0107 *  
okdata2$Edad         -0.08073    0.05345  -1.511   0.1504    
okdata2$SUVR         -0.74099    1.30206  -0.569   0.5772    
okdata2$male         -0.11501    0.77994  -0.147   0.8846    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.562 on 16 degrees of freedom
  (16 observations deleted due to missingness)
Multiple R-squared:  0.4669,	Adjusted R-squared:  0.3003 
F-statistic: 2.803 on 5 and 16 DF,  p-value: 0.05284

> m0 <- lm(okdata2$KD_i ~ okdata2$FPCustom + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR + okdata2$male)
> summary(m0)

Call:
lm(formula = okdata2$KD_i ~ okdata2$FPCustom + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9060 -0.8693  0.0712  1.3062  3.7389 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          67.92106   14.52237   4.677 0.000217 ***
okdata2$FPCustom    -27.68934   41.32332  -0.670 0.511816    
okdata2$Escolaridad   0.31847    0.15604   2.041 0.057094 .  
okdata2$Edad         -0.25900    0.08683  -2.983 0.008356 ** 
okdata2$SUVR          1.23577    2.07456   0.596 0.559237    
okdata2$male          1.14538    1.26034   0.909 0.376169    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.539 on 17 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.4806,	Adjusted R-squared:  0.3278 
F-statistic: 3.146 on 5 and 17 DF,  p-value: 0.03432

> m0 <- lm(okdata2$FName ~ okdata2$FPCustom + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR + okdata2$male)
> summary(m0)

Call:
lm(formula = okdata2$FName ~ okdata2$FPCustom + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.333  -7.319  -1.934   6.399  33.950 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)          143.2128    55.3000   2.590  0.01434 * 
okdata2$FPCustom    -119.3217   158.7761  -0.752  0.45784   
okdata2$Escolaridad    1.3655     0.6463   2.113  0.04251 * 
okdata2$Edad          -1.1730     0.3952  -2.968  0.00564 **
okdata2$SUVR         -11.9182    11.3809  -1.047  0.30285   
okdata2$male          -7.9011     5.6369  -1.402  0.17064   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.39 on 32 degrees of freedom
Multiple R-squared:  0.4612,	Adjusted R-squared:  0.377 
F-statistic: 5.478 on 5 and 32 DF,  p-value: 0.0009417

> m0 <- lm(okdata2$Boston ~ okdata2$FPCustom + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR + okdata2$male)
> summary(m0)

Call:
lm(formula = okdata2$Boston ~ okdata2$FPCustom + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.5752  -1.7398   0.0496   2.4670   7.4385 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          72.9727    15.7856   4.623 6.31e-05 ***
okdata2$FPCustom    -52.9013    45.2722  -1.169  0.25151    
okdata2$Escolaridad   0.5936     0.1842   3.223  0.00298 ** 
okdata2$Edad         -0.2660     0.1129  -2.357  0.02493 *  
okdata2$SUVR          2.6940     3.2501   0.829  0.41350    
okdata2$male         -0.4492     1.6135  -0.278  0.78257    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.386 on 31 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4016,	Adjusted R-squared:  0.305 
F-statistic:  4.16 on 5 and 31 DF,  p-value: 0.005237

> m0 <- lm(okdata2$Boston_Libre ~ okdata2$FPCustom + okdata2$Escolaridad + okdata2$Edad+ okdata2$SUVR + okdata2$male)
> summary(m0)

Call:
lm(formula = okdata2$Boston_Libre ~ okdata2$FPCustom + okdata2$Escolaridad + 
    okdata2$Edad + okdata2$SUVR + okdata2$male)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.6225 -1.9833 -0.2538  1.9214  6.9440 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          76.69651   12.78225   6.000 1.08e-06 ***
okdata2$FPCustom    -57.40968   36.70011  -1.564  0.12759    
okdata2$Escolaridad   0.49401    0.14938   3.307  0.00234 ** 
okdata2$Edad         -0.26234    0.09136  -2.872  0.00719 ** 
okdata2$SUVR          1.92228    2.63062   0.731  0.47026    
okdata2$male         -0.14670    1.30294  -0.113  0.91106    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.558 on 32 degrees of freedom
Multiple R-squared:  0.4376,	Adjusted R-squared:  0.3497 
F-statistic: 4.979 on 5 and 32 DF,  p-value: 0.001748

Resumiendo, En los sujetos con el alelo $\epsilon$-4 presente, hay una ligera asociacion entre la variable LLNaming_total_NP y la fraccion de anisotropia en la Default Mode Network, covariada por edady escolaridad de los sujetos de estudio. Todos los demas efectos que puedan observarse son debido al efecto de asociacion de los resultados de los test con la edad y escolaridad de los sujetos. :-P

Cosas raras

Si hacemos un modelo con todo observamos una cosa curiosa,

> m1 <- lm(okdata$LL_Naming ~ okdata$FPCustom + okdata$SUVR + okdata$Escolaridad + okdata$male +okdata$Edad + okdata$APOE)
> summary(m1)

Call:
lm(formula = okdata$LL_Naming ~ okdata$FPCustom + okdata$SUVR + 
    okdata$Escolaridad + okdata$male + okdata$Edad + okdata$APOE)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.75246  0.02246  0.08062  0.14430  0.30217 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        16.118960   0.507667  31.751   <2e-16 ***
okdata$FPCustom    -2.409148   1.269032  -1.898   0.0596 .  
okdata$SUVR        -0.173373   0.179192  -0.968   0.3349    
okdata$Escolaridad  0.005906   0.005947   0.993   0.3223    
okdata$male        -0.006072   0.056222  -0.108   0.9141    
okdata$Edad        -0.005047   0.003865  -1.306   0.1936    
okdata$APOE        -0.054227   0.045219  -1.199   0.2324    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3216 on 147 degrees of freedom
Multiple R-squared:  0.06926,	Adjusted R-squared:  0.03127 
F-statistic: 1.823 on 6 and 147 DF,  p-value: 0.09836

El modelo es una porqueria pero los residuales indican que hay un estratificacion para la variable LL_Naming,

Y tal y como indica el ajuste,

> m1 <- lm(okdata$LL_Naming ~ okdata$FPCustom)
> summary(m1)

Call:
lm(formula = okdata$LL_Naming ~ okdata$FPCustom)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.88482  0.05739  0.08721  0.12847  0.22546 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      15.5248     0.3728  41.647   <2e-16 ***
okdata$FPCustom  -2.1407     1.2661  -1.691   0.0929 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3248 on 152 degrees of freedom
Multiple R-squared:  0.01846,	Adjusted R-squared:  0.012 
F-statistic: 2.859 on 1 and 152 DF,  p-value: 0.09294

Nota: Estas variables estan muy sesgadas por los valores dados en la clinica, de ahi este comportamiento raro.

Composites

Vamos a intentar el procedimiento utilizando los composites de NP. Las variables de los composites son,

funcioExecutiva_fluencia
funcioExecutiva_velocprocess_IM
funcioExecutiva_atencio
memoria_fnameProf
memoria_fnameNom
memoria_wms
memoria_rbans
gnosia
praxia

Saco las variables que necesito,

> awk -F";" '{print $2,$8,$9,$11,$338,$340,$341,$342,$343,$344,$345,$346,$347,$350,$412}' faceHBI_matriuREF_14-1-19-1.csv | sed 's/ /;/g; s/edat/Edad/; s/subject/Subject/; s/Anyos_Escolaridad_FAC/Escolaridad/; s/Sex_1H_0M/female/; s/Global_v1/SUVR/' > facehbi_data.csv
> awk -F";" {'print $1,$3'} updateAPOE_FACEHBI_051218.csv | sed 's/F//; s/ /,/; s/code_facehbi/Subject/' | sort -n > facehbi_apoe.csv
> sed 's/e2e2\|e2e3/0/; s/e2e4\|e3e3/1/; s/e3e4\|e4e4/2/; s/,/;/' facehbi_apoe.csv > facehbi_apoe_strats.csv
> scp -P 20022 facehbi_apoe_strats.csv detritus.fundacioace.org:facehbi/dti_model/
facehbi_apoe_strats.csv             
> scp -P 20022 facehbi_data.csv detritus.fundacioace.org:facehbi/dti_model/
............
facehbi_[osotolongo@detritus dti_model]$ awk -F";" '{print $1,$4,$6}' facehbi_suvr_dmnfa_customfa_v12.csv | sed 's/ /;/g; s/DMN_FA_v1/DMN/; s/FPCustom_FA_v1/FPCustom/' > facehbi_dti.csv
data.csv

La importo en R,

> fdata <-read.csv("facehbi_data.csv", sep = ";", header=TRUE)
> fapoe <-read.csv("facehbi_apoe_strats.csv", sep = ";", header=TRUE)
> fdti <- read.csv("facehbi_dti.csv", sep=";", header=TRUE)
> okdata <- merge(fdata, fapoe, by = "Subject")
> okdata <- merge(okdata, fdti, by = "Subject")

preparo la lista de variables,

[osotolongo@detritus dti_model]$ cat npvars.names 
funcioExecutiva_fluencia
funcioExecutiva_velocprocess_IM
funcioExecutiva_atencio
memoria_fnameProf
memoria_fnameNom
memoria_wms
memoria_rbans
gnosia
praxia
[osotolongo@detritus dti_model]$ cat nivars.names 
DMN
FPCustom

y hacemos la primera prueba,

> write.csv(okdata, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 2 items

No muy bien hasta ahora,

NP:  funcioExecutiva_velocprocess_IM  NI:  FPCustom
R2:  0.315511475678145  p-value:  0.240278609939977

Estratificamos ahora,

> okdata0 <- okdata[okdata$APOE == "0",]
> write.csv(okdata0, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 2 items
> okdata1 <- okdata[okdata$APOE == "1",]
> write.csv(okdata1, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 2 items
> okdata2 <- okdata[okdata$APOE == "2",]
> write.csv(okdata2, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 2 items

Pero sigue el mismo estilo

APOE 0

NP:  funcioExecutiva_fluencia  NI:  DMN
R2:  0.243867525131838  p-value:  0.362338135897475
NP:  memoria_fnameNom  NI:  DMN
R2:  0.390194941103106  p-value:  0.247812414116152
NP:  memoria_wms  NI:  FPCustom
R2:  0.346181066305969  p-value:  0.61110081152783

APOE 1

NP:  funcioExecutiva_velocprocess_IM  NI:  FPCustom
R2:  0.312707319581829  p-value:  0.159750199653709

APOE 2

NP:  funcioExecutiva_velocprocess_IM  NI:  DMN
R2:  0.438816661720517  p-value:  0.07768386755793
NP:  memoria_fnameNom  NI:  DMN
R2:  0.360648103744986  p-value:  0.680280233588883
NP:  llenguatge_denom_IM  NI:  FPCustom
R2:  0.301108942185496  p-value:  0.168274379413938

y eso es lo mejorcito. Vamos a mirar un poco mejor los APOE 2 (N=38) :-/, que tienen la mejor asociacion.

> m1 <- lm(okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + okdata2$Edad + okdata2$Escolaridad + okdata2$female + okdata2$SUVR)
> summary(m1)
 
Call:
lm(formula = okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + 
    okdata2$Edad + okdata2$Escolaridad + okdata2$female + okdata2$SUVR)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-1.56334 -0.50510  0.08154  0.34458  2.03163 
 
Coefficients:
                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)           1.41530    3.18490   0.444  0.65976   
okdata2$DMN         -14.74565    8.08940  -1.823  0.07768 . 
okdata2$Edad          0.02381    0.01956   1.218  0.23226   
okdata2$Escolaridad  -0.01873    0.02974  -0.630  0.53326   
okdata2$female       -0.45115    0.26240  -1.719  0.09521 . 
okdata2$SUVR          1.68436    0.53322   3.159  0.00345 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.7156 on 32 degrees of freedom
Multiple R-squared:  0.5147,	Adjusted R-squared:  0.4388 
F-statistic: 6.786 on 5 and 32 DF,  p-value: 0.0002047

Voy a limpiar un poco a ver,

 > m1 <- lm(okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + okdata2$female + okdata2$SUVR)
> summary(m1)
 
Call:
lm(formula = okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + 
    okdata2$female + okdata2$SUVR)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-1.72431 -0.45826  0.02665  0.49161  1.97888 
 
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      4.2864     2.3684   1.810 0.079169 .  
okdata2$DMN    -20.8448     6.7607  -3.083 0.004047 ** 
okdata2$female  -0.4298     0.2608  -1.648 0.108592    
okdata2$SUVR     1.9746     0.4927   4.008 0.000317 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.7159 on 34 degrees of freedom
Multiple R-squared:  0.4838,	Adjusted R-squared:  0.4382 
F-statistic: 10.62 on 3 and 34 DF,  p-value: 4.464e-05

un poco mas,

> m1 <- lm(okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + okdata2$SUVR)
> summary(m1)
 
Call:
lm(formula = okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + 
    okdata2$SUVR)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-1.60347 -0.37573 -0.05682  0.41557  2.11191 
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    4.8843     2.3971   2.038 0.049206 *  
okdata2$DMN  -22.6775     6.8300  -3.320 0.002111 ** 
okdata2$SUVR   1.8831     0.5014   3.756 0.000629 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.7333 on 35 degrees of freedom
Multiple R-squared:  0.4426,	Adjusted R-squared:  0.4107 
F-statistic: 13.89 on 2 and 35 DF,  p-value: 3.617e-05

Vaya, no esta tan mal.

A ver si encajo esto de alguna manera,

APOE 0

> m1 <- lm(okdata0$funcioExecutiva_velocprocess_IM ~ okdata0$DMN + okdata0$Escolaridad + okdata0$SUVR)
> summary(m1)

Call:
lm(formula = okdata0$funcioExecutiva_velocprocess_IM ~ okdata0$DMN + 
    okdata0$Escolaridad + okdata0$SUVR)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.88656 -0.46872  0.04536  0.22426  1.55510 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)          2.64148    2.60484   1.014  0.32157   
okdata0$DMN          7.16412    5.08213   1.410  0.17262   
okdata0$Escolaridad -0.11401    0.03127  -3.645  0.00143 **
okdata0$SUVR        -2.79148    2.09112  -1.335  0.19555   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6454 on 22 degrees of freedom
Multiple R-squared:  0.3911,	Adjusted R-squared:  0.3081 
F-statistic:  4.71 on 3 and 22 DF,  p-value: 0.01096

APOE 1

> m1 <- lm(okdata1$funcioExecutiva_velocprocess_IM ~ okdata1$DMN + okdata1$Escolaridad + okdata1$SUVR)
> summary(m1)

Call:
lm(formula = okdata1$funcioExecutiva_velocprocess_IM ~ okdata1$DMN + 
    okdata1$Escolaridad + okdata1$SUVR)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3510 -0.6597 -0.1832  0.3824  4.8138 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)           0.8490     2.0202   0.420 0.675360    
okdata1$DMN          -1.6852     4.8278  -0.349 0.727896    
okdata1$Escolaridad  -0.0851     0.0217  -3.922 0.000176 ***
okdata1$SUVR          0.8174     0.9136   0.895 0.373455    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9785 on 86 degrees of freedom
Multiple R-squared:  0.1598,	Adjusted R-squared:  0.1305 
F-statistic: 5.452 on 3 and 86 DF,  p-value: 0.001767

APOE 2

> m1 <- lm(okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + okdata2$Escolaridad + okdata2$SUVR)
> summary(m1)

Call:
lm(formula = okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$DMN + 
    okdata2$Escolaridad + okdata2$SUVR)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.62197 -0.38893 -0.06813  0.44468  2.06993 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)           4.49262    2.45738   1.828  0.07630 . 
okdata2$DMN         -19.91080    7.67129  -2.595  0.01385 * 
okdata2$Escolaridad  -0.02460    0.03046  -0.808  0.42484   
okdata2$SUVR          1.78264    0.51901   3.435  0.00158 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7369 on 34 degrees of freedom
Multiple R-squared:  0.4531,	Adjusted R-squared:  0.4048 
F-statistic: 9.388 on 3 and 34 DF,  p-value: 0.0001159

Hasta ahora: En los sujetos con el alelo $\epsilon$-4 presente, la variable funcioExecutiva_velocprocess_IM tiene una ligera asociacion con la fraccion de anisotropia en la Default Mode Network y el SUVR de FBB. Esta asociación desaparece cuando no esta presente dicho alelo, siendo la escolaridad del sujeto la variable mas prominente en el caso de los sujetos con el alelo $\epsilon$-2.

Variables de MB

Hay un interes especial en las variables de los tesst “Piramides y Palmeras” y “Kissing and Dancing”. Voy a extraer los valores,

> awk -F";" '{print $2,$8,$9,$11,$286,$287,$288,$289,$412}' faceHBI_matriuREF_14-1-19-1.csv | sed 's/ /;/g; s/edat/Edad/; s/subject/Subject/; s/Anyos_Escolaridad_FAC/Escolaridad/; s/Sex_1H_0M/female/; s/Global_v1/SUVR/; s/Piramides_Y_Palmeras_Palab_FAC/PPp/; s/Piramides_Y_Plameras_Imag_FAC/PPi/; s/Kissing_Dancing_Imagenes_FAC/KDi/; s/Kissing_Dancing_Palabras_FAC/KDp/' > facehbi_data_mb.csv
> scp -P 20022 facehbi_data_mb.csv detritus.fundacioace.org:facehbi/dti_model/
facehbi_data_mb.csv                                                                 100% 5476   801.5KB/s   00:00  

Los cargo, hago un composite con estas variables y miro los modelos.

> fdata <-read.csv("facehbi_data_mb.csv", sep = ";", header=TRUE)
> fapoe <-read.csv("facehbi_apoe_strats.csv", sep = ";", header=TRUE)
> fdti <- read.csv("facehbi_dti.csv", sep=";", header=TRUE)
> okdata <- merge(fdata, fapoe, by = "Subject")
> okdata <- merge(okdata, fdti, by = "Subject")
> okdata$zPPp = (okdata$PPp - mean(okdata$PPp, na.rm = TRUE))/sd(okdata$PPp, na.rm = TRUE)
> okdata$zPPi = (okdata$PPi - mean(okdata$PPi, na.rm = TRUE))/sd(okdata$PPi, na.rm = TRUE)
> okdata$zKDi = (okdata$KDi - mean(okdata$KDi, na.rm = TRUE))/sd(okdata$KDi, na.rm = TRUE)
> okdata$zKDp = (okdata$KDp - mean(okdata$KDp, na.rm = TRUE))/sd(okdata$KDp, na.rm = TRUE)
> np <- data.frame(okdata$zPPp, okdata$zPPi, okdata$zKDp, okdata$zKDi)
> fanp <- fa(np)
> okdata$scop <- fanp$scores
> m1 <- lm(okdata$scop ~ okdata$DMN + okdata$Edad + okdata$Escolaridad + okdata$female + okdata$SUVR)
> summary(m1)

Call:
lm(formula = okdata$scop ~ okdata$DMN + okdata$Edad + okdata$Escolaridad + 
    okdata$female + okdata$SUVR)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7621 -0.2453  0.2214  0.5313  1.0935 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)
(Intercept)         0.79633    2.39102   0.333    0.740
okdata$DMN          3.64776    5.39764   0.676    0.501
okdata$Edad        -0.02282    0.01537  -1.484    0.142
okdata$Escolaridad  0.02556    0.02308   1.108    0.272
okdata$female      -0.19461    0.23710  -0.821    0.415
okdata$SUVR        -0.62693    0.58403  -1.073    0.287

Residual standard error: 0.8891 on 68 degrees of freedom
  (80 observations deleted due to missingness)
Multiple R-squared:  0.1081,	Adjusted R-squared:  0.04254 
F-statistic: 1.649 on 5 and 68 DF,  p-value: 0.1591
> m1 <- lm(okdata$scop ~ okdata$FPCustom + okdata$Edad + okdata$Escolaridad + okdata$female + okdata$SUVR)
> summary(m1)

Call:
lm(formula = okdata$scop ~ okdata$FPCustom + okdata$Edad + okdata$Escolaridad + 
    okdata$female + okdata$SUVR)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8339 -0.2617  0.2257  0.5570  1.1006 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)
(Intercept)         1.33594    2.04679   0.653    0.516
okdata$FPCustom     2.58755    4.99434   0.518    0.606
okdata$Edad        -0.02482    0.01490  -1.666    0.100
okdata$Escolaridad  0.02702    0.02290   1.180    0.242
okdata$female      -0.18170    0.23560  -0.771    0.443
okdata$SUVR        -0.61935    0.58992  -1.050    0.297

Residual standard error: 0.8903 on 68 degrees of freedom
  (80 observations deleted due to missingness)
Multiple R-squared:  0.1057,	Adjusted R-squared:  0.0399 
F-statistic: 1.607 on 5 and 68 DF,  p-value: 0.1701
> okdata0 <- okdata[okdata$APOE == "0",]
> m1 <- lm(okdata0$scop ~ okdata0$FPCustom + okdata0$Edad + okdata0$Escolaridad + okdata0$female + okdata0$SUVR)
> summary(m1)

Call:
lm(formula = okdata0$scop ~ okdata0$FPCustom + okdata0$Edad + 
    okdata0$Escolaridad + okdata0$female + okdata0$SUVR)

Residuals:
      19       20       21       22       23       24       25       26 
 0.85452  0.59106 -0.39878 -1.11895  0.04513  0.21929 -0.25185  0.05957 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)          11.46677   14.83812   0.773    0.520
okdata0$FPCustom    -24.86515   25.51080  -0.975    0.433
okdata0$Edad         -0.09124    0.19281  -0.473    0.683
okdata0$Escolaridad  -0.04443    0.13213  -0.336    0.769
okdata0$female        0.38851    1.52692   0.254    0.823
okdata0$SUVR          2.28677   10.61157   0.215    0.849

Residual standard error: 1.142 on 2 degrees of freedom
  (18 observations deleted due to missingness)
Multiple R-squared:  0.4025,	Adjusted R-squared:  -1.091 
F-statistic: 0.2695 on 5 and 2 DF,  p-value: 0.8972
> okdata1 <- okdata[okdata$APOE == "1",]
> m1 <- lm(okdata1$scop ~ okdata1$FPCustom + okdata1$Edad + okdata1$Escolaridad + okdata1$female + okdata1$SUVR)
> summary(m1)

Call:
lm(formula = okdata1$scop ~ okdata1$FPCustom + okdata1$Edad + 
    okdata1$Escolaridad + okdata1$female + okdata1$SUVR)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7087 -0.3229  0.1671  0.4659  0.7777 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)         -0.709449   2.410371  -0.294    0.770
okdata1$FPCustom     1.379688   4.589883   0.301    0.765
okdata1$Edad        -0.014505   0.013777  -1.053    0.299
okdata1$Escolaridad  0.006632   0.022277   0.298    0.768
okdata1$female      -0.033525   0.229931  -0.146    0.885
okdata1$SUVR         1.141289   1.255969   0.909    0.369

Residual standard error: 0.6507 on 38 degrees of freedom
  (46 observations deleted due to missingness)
Multiple R-squared:  0.04695,	Adjusted R-squared:  -0.07845 
F-statistic: 0.3744 on 5 and 38 DF,  p-value: 0.863
> okdata2 <- okdata[okdata$APOE == "2",]
> m1 <- lm(okdata2$scop ~ okdata2$FPCustom + okdata2$Edad + okdata2$Escolaridad + okdata2$female + okdata2$SUVR)
> summary(m1)

Call:
lm(formula = okdata2$scop ~ okdata2$FPCustom + okdata2$Edad + 
    okdata2$Escolaridad + okdata2$female + okdata2$SUVR)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6076 -0.3908  0.1920  0.7653  1.2746 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)
(Intercept)          2.67538    7.36103   0.363    0.721
okdata2$FPCustom    -2.34987   21.06438  -0.112    0.913
okdata2$Edad        -0.04374    0.04392  -0.996    0.334
okdata2$Escolaridad  0.10479    0.07971   1.315    0.207
okdata2$female      -0.52022    0.64123  -0.811    0.429
okdata2$SUVR        -0.47680    1.07041  -0.445    0.662

Residual standard error: 1.285 on 16 degrees of freedom
  (16 observations deleted due to missingness)
Multiple R-squared:  0.259,	Adjusted R-squared:  0.02737 
F-statistic: 1.118 on 5 and 16 DF,  p-value: 0.3898

Con lo cual me convenzo de que esto no vale para nada.

Todas las redes

No creo que funcione pero por completitud debo hacer el mismo procedimiento para todas las redes que hemos medido, esto es: DMN, FPCustom, LN y SN.

[osotolongo@detritus facehbi]$ awk -F";" '{print $1";"$2}' facehbi_dti_LN.csv > facehbi_fa_LN.csv
[osotolongo@detritus facehbi]$ awk -F";" '{print $1";"$2}' facehbi_dti_SN.csv > facehbi_fa_SN.csv
[osotolongo@detritus facehbi]$ awk -F";" '{print $1";"$2}' facehbi_dti_DMN.csv > facehbi_fa_DMN.csv
[osotolongo@detritus facehbi]$ awk -F";" '{print $1";"$2}' facehbi_dti_FPCustom.csv > facehbi_fa_FPCustom.csv
[osotolongo@detritus facehbi]$ join -t";" -j 1 facehbi_fa_DMN.csv facehbi_fa_LN.csv > tmp.csv; join -t";" -j 1 tmp.csv facehbi_fa_SN.csv > tmp2.csv; join -t";" -j 1 tmp2.csv facehbi_fa_FPCustom.csv > facehbi_fa.csv; rm tmp.csv tmp2.csv
[osotolongo@detritus facehbi]$ head facehbi_fa.csv 
Subject;DMN_FA;LN_FA;SN_FA;FPCustom_FA
001;;;;
002;;;;
003;;;;
004;0.295167;0.297783;0.283152;0.257377
005;0.335169;0.304239;0.326038;0.309156
006;0.306471;0.303277;0.295983;0.285556
007;0.335922;0.323048;0.346638;0.303487
008;0.293183;0.271857;0.279124;0.266333
009;0.343917;0.318176;0.329685;0.319478
[osotolongo@detritus facehbi]$ cp facehbi_fa.csv ~/facehbi/dti_model/

voy a cambiar el script de R para que me de algo mas de info,

library(QuantPsyc)
x<-read.csv("facehbi_dti_np.csv")
Color=c("red","blue")
scan("npvars.names", what = character())->np
scan("nivars.names", what = character())->ni
sink(file = "facehbi_dti_np_models.txt", append = TRUE, type = "output", split = FALSE)
 
for(i in 1:length(np)){
        for(j in 1:length(ni)){
                y.data <- x[c(ni[j], np[i], "female", "Edad", "Escolaridad", "SUVR")]
                y.data <- y.data[complete.cases(y.data),]
                a <- lm( paste ('y.data$', np[i], ' ~ y.data$', ni[j], ' + y.data$SUVR + y.data$female + y.data$Edad + y.data$Escolaridad'))
                writeLines(paste("NP: ", np[i], " NI: ", ni[j]))
                writeLines(paste("R2: ", summary(a)$adj.r.squared, " p-value: ", 1-pf(summary(a)$fstatistic[1], summary(a)$fstatistic[2], summary(a)$fstatistic[3])))
                writeLines(paste("p-value (", ni[j],"): ", summary(a)$coef[2,4], " p-value (SUVR): ", summary(a)$coef[3,4]))
                beta <- lm.beta(a)
                for(k in 1:length(beta)){
                        writeLines(paste(names(beta[k]), ": ", beta[k]))
                }
                writeLines(paste("-------"))
        }
}
sink()

y a probar con los composites de nuevo,

[osotolongo@detritus dti_model]$ cat nivars.names 
DMN_FA
LN_FA
SN_FA
FPCustom_FA
[osotolongo@detritus dti_model]$ cat npvars.names 
funcioExecutiva_fluencia
funcioExecutiva_velocprocess_IM
funcioExecutiva_atencio
memoria_fnameProf
memoria_fnameNom
memoria_wms
memoria_rbans
gnosia
praxia
llenguatge_denom_IM

Empiezo con el global,

> fdti <- read.csv("facehbi_fa.csv", sep=";", header=TRUE)
> fdata <-read.csv("facehbi_data.csv", sep = ";", header=TRUE)
> fapoe <-read.csv("facehbi_apoe_strats.csv", sep = ";", header=TRUE)
> okdata <- merge(fdata, fapoe, by = "Subject")
> okdata <- merge(okdata, fdti, by = "Subject")
> write.csv(okdata, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 4 items

movemos los resultados,

[osotolongo@detritus dti_model]$ mv facehbi_dti_np_models.txt facehbi_dti_np_models_all.txt

y ahora a estratificar,

> okdata0 <- okdata[okdata$APOE == "0",]
> write.csv(okdata0, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 4 items
[osotolongo@detritus dti_model]$ mv facehbi_dti_np_models.txt facehbi_dti_np_models_0.txt
> okdata1 <- okdata[okdata$APOE == "1",]
> write.csv(okdata1, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 4 items
[osotolongo@detritus dti_model]$ mv facehbi_dti_np_models.txt facehbi_dti_np_models_1.txt
> okdata2 <- okdata[okdata$APOE == "2",]
> write.csv(okdata2, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 10 items
Read 4 items
[osotolongo@detritus dti_model]$ mv facehbi_dti_np_models.txt facehbi_dti_np_models_2.txt

y lo voy a guardar, por si acaso,

[osotolongo@detritus dti_model]$ tar czvf facehbi_dti_np_models.tgz facehbi_dti_np_models_*
facehbi_dti_np_models_0.txt
facehbi_dti_np_models_1.txt
facehbi_dti_np_models_2.txt
facehbi_dti_np_models_all.txt

Nota: Estos hay que revisarlos despacio pues puede haber alguna asosiacion con el SUVR que hayamos pasado por alto.

Ahora el otro composite,

> fdata <-read.csv("facehbi_data_mb.csv", sep = ";", header=TRUE)
> okdata <- merge(okdata, fdata, by = "Subject")
> okdata$zPPp = (okdata$PPp - mean(okdata$PPp, na.rm = TRUE))/sd(okdata$PPp, na.rm = TRUE)
> okdata$zPPi = (okdata$PPi - mean(okdata$PPi, na.rm = TRUE))/sd(okdata$PPi, na.rm = TRUE)
> okdata$zKDi = (okdata$KDi - mean(okdata$KDi, na.rm = TRUE))/sd(okdata$KDi, na.rm = TRUE)
> okdata$zKDp = (okdata$KDp - mean(okdata$KDp, na.rm = TRUE))/sd(okdata$KDp, na.rm = TRUE)
> mb <- data.frame(okdata$zPPp, okdata$zPPi, okdata$zKDp, okdata$zKDi)
> okdata$cs <- mbsc$scores
[osotolongo@detritus dti_model]$ cat npvars.names
cs
> write.csv(okdata, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 1 item
Read 4 items

y lo hago estratificado tambien, (moviendo los outputs como antes)

> okdata0 <- okdata[okdata$APOE == "0",]
> write.csv(okdata0, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 1 item
Read 4 items
> okdata1 <- okdata[okdata$APOE == "1",]
> write.csv(okdata1, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 1 item
Read 4 items
> okdata2 <- okdata[okdata$APOE == "2",]
> write.csv(okdata2, file="facehbi_dti_np.csv")
> source("get_lms.r")
Read 1 item
Read 4 items

y aqui si no hay nada de nada. Voy a hacerme un script para sacar cuando R2 es mayor que 0.3 por poner un numero,

checkr2.pl
#!/usr/bin/perl
 
use strict;
use warnings;
use Data::Dump qw(dump);
 
my $ifile = "facehbi_dti_np_models_all.txt";
my $thresh = 0.3;
my %model;
open IDF, "<$ifile" or die "No such file\n";
while (<IDF>){
        if (/-------/ && $model{"r2"}>0.3) {
                print $model{"ni_var"}.",".$model{np_var}."\nr2 = ".$model{"r2"}.", p-value = ".$model{"pvalue"}."\npv_".$model{"ni_var"}." = ".$model{"pv_ni"}.", pv_SUVR = ".$model{pv_suvr}." \n";
        };
        if (/^NP:.*/) {($model{"np_var"}, $model{"ni_var"}) = /^NP:\s+(\w+)\s+NI:\s+(\w+)\s*$/};
        if (/^R2:.*/) {($model{"r2"}, $model{"pvalue"}) = /^R2:\s+(\S+)\s+p-value:\s+(\S+)/};
        if (/^p-value.*/) {($model{"pv_ni"}, $model{"pv_suvr"}) = /p-value.*:\s+(\S+)\s+p-value.*:\s+(\S+)/}
}
close IDF;

Para todos,

[osotolongo@detritus dti_model]$ ./checkr2.pl 
Analizing facehbi_dti_np_models_all.txt ...

DMN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.302226155665606, p-value = 1.59872115546023e-14
pv_DMN_FA = 0.883576950793585, pv_SUVR = 0.24424229051191 

LN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.302579788141454, p-value = 1.53210777398272e-14
pv_LN_FA = 0.731754149108559, pv_SUVR = 0.239811378742586 

SN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.302699583800708, p-value = 1.50990331349021e-14
pv_SN_FA = 0.698477685239432, pv_SUVR = 0.241344996365449 

FPCustom_FA, funcioExecutiva_velocprocess_IM
r2 = 0.303857285316447, p-value = 1.28785870856518e-14
pv_FPCustom_FA = 0.495322950618766, pv_SUVR = 0.229845691710707 

APOE 0

[osotolongo@detritus dti_model]$ ./checkr2.pl 
Analizing facehbi_dti_np_models_0.txt ...

DMN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.379552223636583, p-value = 0.0057104958050207
pv_DMN_FA = 0.639790744923555, pv_SUVR = 0.0157173219667127 

LN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.373746758451691, p-value = 0.00627788677009611
pv_LN_FA = 0.922726301726418, pv_SUVR = 0.0179719566898906 

SN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.373504904382816, p-value = 0.00630256049240607
pv_SN_FA = 0.978604565669928, pv_SUVR = 0.0176943724392262 

FPCustom_FA, funcioExecutiva_velocprocess_IM
r2 = 0.376679557258316, p-value = 0.00598535915763587
pv_FPCustom_FA = 0.734462972707933, pv_SUVR = 0.0164897973831555 

DMN_FA, memoria_wms
r2 = 0.384900678414176, p-value = 0.00522814187347165
pv_DMN_FA = 0.971755525011924, pv_SUVR = 0.126202366929838 

LN_FA, memoria_wms
r2 = 0.476189288707085, p-value = 0.000987730906959694
pv_LN_FA = 0.0571602758952463, pv_SUVR = 0.0786860000743259 

SN_FA, memoria_wms
r2 = 0.4105573407532, p-value = 0.00337775223148928
pv_SN_FA = 0.327134537994976, pv_SUVR = 0.110122034640346 

FPCustom_FA, memoria_wms
r2 = 0.38938771125582, p-value = 0.00485148460500173
pv_FPCustom_FA = 0.683663237630694, pv_SUVR = 0.141496294121895 

APOE 1 –> nada

APOE 2

[osotolongo@detritus dti_model]$ ./checkr2.pl 
Analizing facehbi_dti_np_models_2.txt ...

DMN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.485798156550142, p-value = 2.68765592381648e-06
pv_DMN_FA = 0.00443564461139048, pv_SUVR = 0.00250803493024018 

LN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.391142446695136, p-value = 7.40330891934038e-05
pv_LN_FA = 0.298247246087403, pv_SUVR = 0.00913580342963226 

SN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.473679765360014, p-value = 4.26313995949279e-06
pv_SN_FA = 0.0075840881661726, pv_SUVR = 0.00315334075632377 

FPCustom_FA, funcioExecutiva_velocprocess_IM
r2 = 0.420481875101633, p-value = 2.83138960435192e-05
pv_FPCustom_FA = 0.0767718841930776, pv_SUVR = 0.00977255081753087 

DMN_FA, memoria_fnameProf
r2 = 0.327588088792215, p-value = 0.000497956057757376
pv_DMN_FA = 0.773974302319369, pv_SUVR = 0.788744944522824 

LN_FA, memoria_fnameProf
r2 = 0.382914152316493, p-value = 9.60047895888216e-05
pv_LN_FA = 0.0561917139565343, pv_SUVR = 0.624162892962772 

SN_FA, memoria_fnameProf
r2 = 0.327676768593875, p-value = 0.000496710477748463
pv_SN_FA = 0.76680425794751, pv_SUVR = 0.789799021308529 

FPCustom_FA, memoria_fnameProf
r2 = 0.331099185559454, p-value = 0.000450805039077462
pv_FPCustom_FA = 0.584035234528739, pv_SUVR = 0.810357038625016 

DMN_FA, memoria_fnameNom
r2 = 0.328367049483929, p-value = 0.000487113972439834
pv_DMN_FA = 0.898296424312723, pv_SUVR = 0.0560073342550184 

LN_FA, memoria_fnameNom
r2 = 0.330867944990235, p-value = 0.000453777260227328
pv_LN_FA = 0.679077371093828, pv_SUVR = 0.0600344325872235 

SN_FA, memoria_fnameNom
r2 = 0.331778400581458, p-value = 0.000442179712075386
pv_SN_FA = 0.633252108652689, pv_SUVR = 0.0601087607629487 

FPCustom_FA, memoria_fnameNom
r2 = 0.334618953643581, p-value = 0.000407752047013266
pv_FPCustom_FA = 0.524784142484717, pv_SUVR = 0.0519653979539324 

DMN_FA, memoria_wms
r2 = 0.30272247894919, p-value = 0.000989109715828107
pv_DMN_FA = 0.228059068988977, pv_SUVR = 0.472409952167858 

Voy a mirar un poco,

> m0 <- lm(okdata0$funcioExecutiva_velocprocess_IM ~ okdata0$SUVR + okdata0$Edad)
> summary(m0)
 
Call:
lm(formula = okdata0$funcioExecutiva_velocprocess_IM ~ okdata0$SUVR + 
    okdata0$Edad)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-1.3337 -0.5111 -0.1135  0.3141  2.5369 
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept)   2.79206    3.21204   0.869  0.39267   
okdata0$SUVR -6.70012    2.41598  -2.773  0.01013 * 
okdata0$Edad  0.07533    0.02137   3.526  0.00159 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.9002 on 26 degrees of freedom
Multiple R-squared:  0.4416,	Adjusted R-squared:  0.3987 
F-statistic: 10.28 on 2 and 26 DF,  p-value: 0.0005129
 
> m1 <- lm(okdata1$funcioExecutiva_velocprocess_IM ~ okdata1$SUVR + okdata1$Edad)
> summary(m1)
 
Call:
lm(formula = okdata1$funcioExecutiva_velocprocess_IM ~ okdata1$SUVR + 
    okdata1$Edad)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-1.2138 -0.5730 -0.1959  0.2146  4.7844 
 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.200673   1.052358  -3.992 0.000114 ***
okdata1$SUVR  0.002109   0.780712   0.003 0.997849    
okdata1$Edad  0.063894   0.011991   5.328 4.82e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.8937 on 118 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.2052,	Adjusted R-squared:  0.1918 
F-statistic: 15.24 on 2 and 118 DF,  p-value: 1.3e-06
 
> m2 <- lm(okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$SUVR + okdata2$Edad)
> summary(m2)
 
Call:
lm(formula = okdata2$funcioExecutiva_velocprocess_IM ~ okdata2$SUVR + 
    okdata2$Edad)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-1.74708 -0.28264 -0.07348  0.26384  2.76441 
 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.49416    1.14875  -3.912   0.0003 ***
okdata2$SUVR  1.34940    0.54032   2.497   0.0161 *  
okdata2$Edad  0.04144    0.01738   2.384   0.0213 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.8034 on 46 degrees of freedom
Multiple R-squared:  0.2749,	Adjusted R-squared:  0.2434 
F-statistic: 8.719 on 2 and 46 DF,  p-value: 0.0006158

Riesgo - No riesgo

Vamosaplantear el problema de manera distinta. Supongamos que la contribucion del APOE depende solo de la presencia del alelo $\epsilon$-4 y clasifiquemos los sujetos segun esto, en con riesgo o sin riesgo.

> okdata$Risk <- ifelse (okdata$APOE==2 , 1, 0)

Pero ahora voy a hacer una cosa un poco mas complicada,

get_lms2.r
library(QuantPsyc)
x<-read.csv("facehbi_dti_np.csv")
Color=c("red","blue")
scan("npvars.names", what = character())->np
scan("nivars.names", what = character())->ni
sink(file = "facehbi_dti_np_models.txt", append = TRUE, type = "output", split = FALSE)
 
for(i in 1:length(np)){
        for(j in 1:length(ni)){
                y.data <- x[c(ni[j], np[i], "female", "Edad", "Escolaridad", "SUVR", "Risk")]
                y.data <- y.data[complete.cases(y.data),]
                a <- lm( paste ('y.data$', np[i], ' ~ y.data$', ni[j], ' + y.data$SUVR +y.data$Risk + y.data$female + y.data$Edad + y.data$Escolaridad + ', 'y.data$', ni[j], '*y.data$Risk'))
                writeLines(paste("NP: ", np[i], " NI: ", ni[j]))
                writeLines(paste("R2: ", summary(a)$adj.r.squared, " p-value: ", 1-pf(summary(a)$fstatistic[1], summary(a)$fstatistic[2], summary(a)$fstatistic[3])))
                writeLines(paste("p-value (", ni[j],"): ", summary(a)$coef[2,4], " p-value (SUVR): ", summary(a)$coef[3,4]))
                beta <- lm.beta(a)
                for(k in 1:length(beta)){
                        writeLines(paste(names(beta[k]), ": ", beta[k]))
                }
                writeLines(paste("-------"))
        }
}
sink()

Asi que pruebo con el global,

> write.csv(okdata, file="facehbi_dti_np.csv")
> source("get_lms2.r")

y luego,

[osotolongo@detritus dti_model]$ ./checkr2.pl 
Analizing facehbi_dti_np_models.txt ...
 
DMN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.312799089824939, p-value = 2.93098878501041e-14
pv_DMN_FA = 0.476462254461098, pv_SUVR = 0.179281559330912 
 
SN_FA, funcioExecutiva_velocprocess_IM
r2 = 0.311504329049638, p-value = 3.47499806707674e-14
pv_SN_FA = 0.551975051522526, pv_SUVR = 0.181998130409768 
 
FPCustom_FA, funcioExecutiva_velocprocess_IM
r2 = 0.311678871879767, p-value = 3.39728245535298e-14
pv_FPCustom_FA = 0.151421806156447, pv_SUVR = 0.201473541858229 

puaf, a ver,

> m <- lm(okdata$funcioExecutiva_velocprocess_IM ~ okdata$SUVR + okdata$Edad + okdata$Escolaridad + okdata$female + okdata$DMN_FA*okdata$Risk)
> summary(m)
 
Call:
lm(formula = okdata$funcioExecutiva_velocprocess_IM ~ okdata$SUVR + 
    okdata$Edad + okdata$Escolaridad + okdata$female + okdata$DMN_FA * 
    okdata$Risk)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-1.4094 -0.5672 -0.1264  0.3442  4.3680 
 
Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                -4.12971    1.33068  -3.103  0.00221 ** 
okdata$SUVR                 0.58201    0.43176   1.348  0.17928    
okdata$Edad                 0.05310    0.00865   6.139 4.82e-09 ***
okdata$Escolaridad         -0.04301    0.01362  -3.157  0.00186 ** 
okdata$female              -0.38688    0.13066  -2.961  0.00346 ** 
okdata$DMN_FA               2.24716    3.14978   0.713  0.47646    
okdata$Risk                 5.03080    2.38260   2.111  0.03605 *  
okdata$DMN_FA:okdata$Risk -15.70430    7.30118  -2.151  0.03276 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.8286 on 188 degrees of freedom
  (4 observations deleted due to missingness)
Multiple R-squared:  0.3375,	Adjusted R-squared:  0.3128 
F-statistic: 13.68 on 7 and 188 DF,  p-value: 2.936e-14

No, gracias. :-\

neuroimagen/altdti.txt · Last modified: 2020/08/04 10:58 (external edit)