ModelsPPS

#Importing

library(haven)
pps_data_dti <- read_sav("paraps_data_dti_fasci_tagliati.sav")
pps <- as.data.frame(pps_data_dti)
pps$group <- as.factor(pps$group)
#rinomino UPPERLIMB in UL e LOWER in LL
names <- names(pps)
names<-gsub("UPPERLIMB", "UL",names)
names<-gsub("LOWERLIMB", "LL",names)
names

  [1] "subj"          "group"         "age"           "sex"          
  [5] "LL_CP_LH_FA"   "LL_CP_LH_MD"   "LL_CP_LH_AD"   "LL_CP_LH_RD"  
  [9] "LL_CP_RH_FA"   "LL_CP_RH_MD"   "LL_CP_RH_AD"   "LL_CP_RH_RD"  
 [13] "UL_CP_LH_FA"   "UL_CP_LH_MD"   "UL_CP_LH_AD"   "UL_CP_LH_RD"  
 [17] "UL_CP_RH_FA"   "UL_CP_RH_MD"   "UL_CP_RH_AD"   "UL_CP_RH_RD"  
 [21] "SMA_CP_LH_FA"  "SMA_CP_LH_MD"  "SMA_CP_LH_AD"  "SMA_CP_LH_RD" 
 [25] "SMA_CP_RH_FA"  "SMA_CP_RH_MD"  "SMA_CP_RH_AD"  "SMA_CP_RH_RD" 
 [29] "LL_VLA_LH_FA"  "LL_VLA_LH_MD"  "LL_VLA_LH_AD"  "LL_VLA_LH_RD" 
 [33] "LL_VLA_RH_FA"  "LL_VLA_RH_MD"  "LL_VLA_RH_AD"  "LL_VLA_RH_RD" 
 [37] "UL_VLA_LH_FA"  "UL_VLA_LH_MD"  "UL_VLA_LH_AD"  "UL_VLA_LH_RD" 
 [41] "UL_VLA_RH_FA"  "UL_VLA_RH_MD"  "UL_VLA_RH_AD"  "UL_VLA_RH_RD" 
 [45] "SMA_VLA_LH_FA" "SMA_VLA_LH_MD" "SMA_VLA_LH_AD" "SMA_VLA_LH_RD"
 [49] "SMA_VLA_RH_FA" "SMA_VLA_RH_MD" "SMA_VLA_RH_AD" "SMA_VLA_RH_RD"
 [53] "PMA_VLA_LH_FA" "PMA_VLA_LH_MD" "PMA_VLA_LH_AD" "PMA_VLA_LH_RD"
 [57] "PMA_VLA_RH_FA" "PMA_VLA_RH_MD" "PMA_VLA_RH_AD" "PMA_VLA_RH_RD"
 [61] "LL_VLP_LH_FA"  "LL_VLP_LH_MD"  "LL_VLP_LH_AD"  "LL_VLP_LH_RD" 
 [65] "LL_VLP_RH_FA"  "LL_VLP_RH_MD"  "LL_VLP_RH_AD"  "LL_VLP_RH_RD" 
 [69] "UL_VLP_LH_FA"  "UL_VLP_LH_MD"  "UL_VLP_LH_AD"  "UL_VLP_LH_RD" 
 [73] "UL_VLP_RH_FA"  "UL_VLP_RH_MD"  "UL_VLP_RH_AD"  "UL_VLP_RH_RD" 
 [77] "SMA_VLP_LH_FA" "SMA_VLP_LH_MD" "SMA_VLP_LH_AD" "SMA_VLP_LH_RD"
 [81] "SMA_VLP_RH_FA" "SMA_VLP_RH_MD" "SMA_VLP_RH_AD" "SMA_VLP_RH_RD"
 [85] "PMA_VLP_LH_FA" "PMA_VLP_LH_MD" "PMA_VLP_LH_AD" "PMA_VLP_LH_RD"
 [89] "PMA_VLP_RH_FA" "PMA_VLP_RH_MD" "PMA_VLP_RH_AD" "PMA_VLP_RH_RD"
 [93] "V1_PUL_LH_FA"  "V1_PUL_LH_MD"  "V1_PUL_LH_AD"  "V1_PUL_LH_RD" 
 [97] "V1_PUL_RH_FA"  "V1_PUL_RH_MD"  "V1_PUL_RH_AD"  "V1_PUL_RH_RD" 
[101] "filter_$"     

names(pps) <- names

library(nlme);library("ez"); library("ggplot2"); library("nlme");library("pastecs"); library("reshape");library("WRS2")

CP

FA

pps_CP_covs_FA <- pps[,c("subj","group","sex", "age","LL_CP_LH_FA", "LL_CP_RH_FA", "UL_CP_LH_FA", "UL_CP_RH_FA")]
CP_covs_FA_l <- melt(pps_CP_covs_FA, id = c("subj", "group", "sex", "age"), measured = c("LL_CP_LH_FA", "LL_CP_RH_FA", "UL_CP_LH_FA", "UL_CP_RH_FA"))
names(CP_covs_FA_l)[6] <- "FA"
CP_covs_FA_l$hemi <- gl(2, 80, labels = c("L-hemi", "R-hemi"))
CP_covs_FA_l$limb <- gl(2, 160, labels = c("LOWER", "UPPER"))

CP_covs_FA_l <- CP_covs_FA_l[-which(is.na(CP_covs_FA_l$FA)),]

#model
mod.CP_grp_limb_covs_FA <- lme(FA ~ -1+group*limb + hemi + sex + age,  random = ~1| subj, data = CP_covs_FA_l, method = "ML")
summary(mod.CP_grp_limb_covs_FA)

Linear mixed-effects model fit by maximum likelihood
  Data: CP_covs_FA_l 
        AIC       BIC   logLik
  -1269.679 -1235.793 643.8397

Random effects:
 Formula: ~1 | subj
        (Intercept)   Residual
StdDev:  0.03982907 0.02341806

Fixed effects:  FA ~ -1 + group * limb + hemi + sex + age 
                      Value   Std.Error  DF   t-value p-value
group0            0.4861571 0.020376947  76 23.858192  0.0000
group1            0.4517929 0.020294674  76 22.261649  0.0000
limbUPPER        -0.0061432 0.003744027 237 -1.640807  0.1022
hemiR-hemi       -0.0088420 0.002652867 237 -3.332985  0.0010
sex               0.0081552 0.010130440  76  0.805021  0.4233
age              -0.0005836 0.000332841  76 -1.753437  0.0836
group1:limbUPPER  0.0084232 0.005305734 237  1.587571  0.1137
 Correlation: 
                 group0 group1 lUPPER hmR-hm sex    age   
group1            0.885                                   
limbUPPER        -0.092  0.000                            
hemiR-hemi       -0.064 -0.064  0.000                     
sex              -0.457 -0.458  0.000  0.001              
age              -0.878 -0.877  0.000 -0.001  0.147       
group1:limbUPPER  0.064 -0.067 -0.706 -0.004 -0.001  0.001

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.07013059 -0.52295641 -0.03112463  0.59214241  2.40286648 

Number of Observations: 319
Number of Groups: 80 

Residuals

#resid(intermod)
hist(resid(mod.CP_grp_limb_covs_FA))

shapiro.test(resid(mod.CP_grp_limb_covs_FA))

    Shapiro-Wilk normality test

data:  resid(mod.CP_grp_limb_covs_FA)
W = 0.99308, p-value = 0.1487

MultComp correction

library(multcomp)

Caricamento del pacchetto richiesto: mvtnorm

Caricamento del pacchetto richiesto: survival

Caricamento del pacchetto richiesto: TH.data

Caricamento del pacchetto richiesto: MASS

Caricamento pacchetto: 'TH.data'

Il seguente oggetto è mascherato da 'package:MASS':

    geyser

# library()
# contr <- rbind( "UPLC - UPLP" = c(1, 0, 0, 0,-1, 0, 0, 0),
#                 "LOLC - LOLP" = c(0, 1, 0, 0, 0,-1, 0, 0),
#                 "UPRC - UPRP" = c(0, 0, 1, 0, 0, 0,-1, 0),
#                 "LORC - LORP" = c(0, 0, 0, 1, 0, 0, 0,-1))
# #ask : inverse contrasts?
# 
# mcontr.CP_grp_limb_covs_FA=glht(mod.CP_grp_limb_covs_FA, linfct = contr)
# summary(mcontr.CP_grp_limb_covs_FA)

MD

pps_CP_covs_MD <- pps[,c("subj","group","sex", "age","LL_CP_LH_MD", "LL_CP_RH_MD", "UL_CP_LH_MD", "UL_CP_RH_MD")]
CP_covs_MD_l <- melt(pps_CP_covs_MD, id = c("subj", "group", "sex", "age"), measured = c("LL_CP_LH_MD", "LL_CP_RH_MD", "UL_CP_LH_MD", "UL_CP_RH_MD"))
names(CP_covs_MD_l)[6] <- "MD"
CP_covs_MD_l$hemi <- gl(2, 80, labels = c("L-hemi", "R-hemi"))
CP_covs_MD_l$limb <- gl(2, 160, labels = c("LOWER", "UPPER"))

CP_covs_MD_l <- CP_covs_MD_l[-which(is.na(CP_covs_MD_l$MD)),]

#model
mod.CP_grp_limb_covs_MD <- lme(MD ~ -1+group*limb + hemi + sex + age,  random = ~1| subj, data = CP_covs_MD_l, method = "ML")
summary(mod.CP_grp_limb_covs_MD)

Linear mixed-effects model fit by maximum likelihood
  Data: CP_covs_MD_l 
        AIC       BIC   logLik
  -5441.646 -5407.759 2729.823

Random effects:
 Formula: ~1 | subj
         (Intercept)     Residual
StdDev: 3.168489e-05 3.966618e-05

Fixed effects:  MD ~ -1 + group * limb + hemi + sex + age 
                         Value    Std.Error  DF  t-value p-value
group0            0.0007194811 1.864137e-05  76 38.59593  0.0000
group1            0.0007245210 1.857946e-05  76 38.99581  0.0000
limbUPPER        -0.0000042500 6.341740e-06 237 -0.67016  0.5034
hemiR-hemi        0.0000132789 4.492951e-06 237  2.95549  0.0034
sex              -0.0000007749 9.122605e-06  76 -0.08494  0.9325
age               0.0000001907 2.998070e-07  76  0.63604  0.5267
group1:limbUPPER -0.0000069578 8.985902e-06 237 -0.77430  0.4395
 Correlation: 
                 group0 group1 lUPPER hmR-hm sex    age   
group1            0.867                                   
limbUPPER        -0.170  0.000                            
hemiR-hemi       -0.119 -0.118  0.000                     
sex              -0.450 -0.450  0.000  0.001              
age              -0.864 -0.863  0.000 -0.002  0.146       
group1:limbUPPER  0.119 -0.123 -0.706 -0.004 -0.001  0.002

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.7666346 -0.5816299 -0.1236227  0.4144602  3.7630666 

Number of Observations: 319
Number of Groups: 80 

VLA

FA

pps_VLA_covs_FA <- pps[,c("subj","group","sex", "age","LL_VLA_LH_FA", "LL_VLA_RH_FA", "UL_VLA_LH_FA", "UL_VLA_RH_FA")]
VLA_covs_FA_l <- melt(pps_VLA_covs_FA, id = c("subj", "group", "sex", "age"), measured = c("LL_VLA_LH_FA", "LL_VLA_RH_FA", "UL_VLA_LH_FA", "UL_VLA_RH_FA"))
names(VLA_covs_FA_l)[6] <- "FA"
VLA_covs_FA_l$hemi <- gl(2, 80, labels = c("L-hemi", "R-hemi"))
VLA_covs_FA_l$limb <- gl(2, 160, labels = c("LOWER", "UPPER"))

VLA_covs_FA_l <- VLA_covs_FA_l[-which(is.na(VLA_covs_FA_l$FA)),]

mod.VLA_grp_limb_covs_FA <- lme(FA ~ -1+group*limb + hemi + sex + age,  random = ~1| subj, data = VLA_covs_FA_l, method = "ML")
summary(mod.VLA_grp_limb_covs_FA)

Linear mixed-effects model fit by maximum likelihood
  Data: VLA_covs_FA_l 
        AIC       BIC   logLik
  -892.7519 -859.5997 455.3759

Random effects:
 Formula: ~1 | subj
        (Intercept)   Residual
StdDev:  0.04096649 0.04198365

Fixed effects:  FA ~ -1 + group * limb + hemi + sex + age 
                      Value   Std.Error  DF   t-value p-value
group0            0.5138993 0.023031510  76 22.312877  0.0000
group1            0.4807230 0.022927292  76 20.967282  0.0000
limbUPPER         0.0105510 0.007229903 212  1.459362  0.1459
hemiR-hemi       -0.0008274 0.005001953 212 -0.165423  0.8688
sex              -0.0096406 0.011372544  76 -0.847712  0.3993
age              -0.0007469 0.000371767  76 -2.008935  0.0481
group1:limbUPPER  0.0005403 0.010077966 212  0.053616  0.9573
 Correlation: 
                 group0 group1 lUPPER hmR-hm sex    age   
group1            0.867                                   
limbUPPER        -0.158  0.011                            
hemiR-hemi       -0.095 -0.090  0.035                     
sex              -0.449 -0.449 -0.012 -0.015              
age              -0.864 -0.867 -0.014 -0.015  0.145       
group1:limbUPPER  0.109 -0.123 -0.718 -0.028  0.008  0.016

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-9.08415811 -0.30287509 -0.03885641  0.34021435  3.89966938 

Number of Observations: 294
Number of Groups: 80 

MD

pps_VLA_covs_MD <- pps[,c("subj","group","sex", "age","LL_VLA_LH_MD", "LL_VLA_RH_MD", "UL_VLA_LH_MD", "UL_VLA_RH_MD")]
VLA_covs_MD_l <- melt(pps_VLA_covs_MD, id = c("subj", "group", "sex", "age"), measured = c("LL_VLA_LH_MD", "LL_VLA_RH_MD", "UL_VLA_LH_MD", "UL_VLA_RH_MD"))
names(VLA_covs_MD_l)[6] <- "MD"
VLA_covs_MD_l$hemi <- gl(2, 80, labels = c("L-hemi", "R-hemi"))
VLA_covs_MD_l$limb <- gl(2, 160, labels = c("LOWER", "UPPER"))

VLA_covs_MD_l <- VLA_covs_MD_l[-which(is.na(VLA_covs_MD_l$MD)),]

mod.VLA_grp_limb_covs_MD <- lme(MD ~ -1+group*limb + hemi + sex + age,  random = ~1| subj, data = VLA_covs_MD_l, method = "ML")
summary(mod.VLA_grp_limb_covs_MD)

Linear mixed-effects model fit by maximum likelihood
  Data: VLA_covs_MD_l 
        AIC       BIC   logLik
  -4964.465 -4931.344 2491.233

Random effects:
 Formula: ~1 | subj
         (Intercept)     Residual
StdDev: 3.235102e-05 4.201712e-05

Fixed effects:  MD ~ -1 + group * limb + hemi + sex + age 
                         Value    Std.Error  DF  t-value p-value
group0            0.0007131003 1.948724e-05  76 36.59319  0.0000
group1            0.0007255123 1.942606e-05  76 37.34738  0.0000
limbUPPER        -0.0000204152 7.265737e-06 211 -2.80979  0.0054
hemiR-hemi        0.0000087625 5.005117e-06 211  1.75071  0.0814
sex               0.0000136802 9.573548e-06  76  1.42895  0.1571
age               0.0000000927 3.125670e-07  76  0.29655  0.7676
group1:limbUPPER -0.0000032496 1.009818e-05 211 -0.32180  0.7479
 Correlation: 
                 group0 group1 lUPPER hmR-hm sex    age   
group1            0.857                                   
limbUPPER        -0.184  0.019                            
hemiR-hemi       -0.113 -0.109  0.028                     
sex              -0.444 -0.445 -0.018 -0.016              
age              -0.856 -0.861 -0.021 -0.015  0.146       
group1:limbUPPER  0.127 -0.148 -0.720 -0.022  0.012  0.022

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.3071972 -0.5313491 -0.1110331  0.4308241  5.5740491 

Number of Observations: 293
Number of Groups: 80 

VLP

FA

pps_VLP_covs_FA <- pps[,c("subj","group","sex", "age","LL_VLP_LH_FA", "LL_VLP_RH_FA", "UL_VLP_LH_FA", "UL_VLP_RH_FA")]
VLP_covs_FA_l <- melt(pps_VLP_covs_FA, id = c("subj", "group", "sex", "age"), measured = c("LL_VLP_LH_FA", "LL_VLP_RH_FA", "UL_VLP_LH_FA", "UL_VLP_RH_FA"))
names(VLP_covs_FA_l)[6] <- "FA"
VLP_covs_FA_l$hemi <- gl(2, 80, labels = c("L-hemi", "R-hemi"))
VLP_covs_FA_l$limb <- gl(2, 160, labels = c("LOWER", "UPPER"))

VLP_covs_FA_l <- VLP_covs_FA_l[-which(is.na(VLP_covs_FA_l$FA)),]

#model
mod.VLP_grp_limb_covs_FA <- lme(FA ~ -1+group*limb + hemi + sex + age,  random = ~1| subj, data = VLP_covs_FA_l, method = "ML")
summary(mod.VLP_grp_limb_covs_FA)

Linear mixed-effects model fit by maximum likelihood
  Data: VLP_covs_FA_l 
        AIC       BIC   logLik
  -1225.168 -1191.337 621.5838

Random effects:
 Formula: ~1 | subj
        (Intercept)   Residual
StdDev:  0.04388103 0.02447551

Fixed effects:  FA ~ -1 + group * limb + hemi + sex + age 
                      Value   Std.Error  DF   t-value p-value
group0            0.4856864 0.022352778  76 21.728235  0.0000
group1            0.4496527 0.022280961  76 20.181029  0.0000
limbUPPER         0.0035115 0.003913366 235  0.897316  0.3705
hemiR-hemi       -0.0091757 0.002781611 235 -3.298695  0.0011
sex               0.0091680 0.011121514  76  0.824352  0.4123
age              -0.0005617 0.000365362  76 -1.537359  0.1284
group1:limbUPPER  0.0093684 0.005579272 235  1.679148  0.0945
 Correlation: 
                 group0 group1 lUPPER hmR-hm sex    age   
group1            0.885                                   
limbUPPER        -0.088  0.000                            
hemiR-hemi       -0.062 -0.061  0.000                     
sex              -0.457 -0.458  0.000 -0.002              
age              -0.879 -0.877  0.000  0.000  0.146       
group1:limbUPPER  0.059 -0.068 -0.701 -0.004  0.000  0.003

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.25963916 -0.49104745 -0.03899312  0.52735891  2.64428583 

Number of Observations: 317
Number of Groups: 80 

MD

pps_VLP_covs_MD <- pps[,c("subj","group","sex", "age","LL_VLP_LH_MD", "LL_VLP_RH_MD", "UL_VLP_LH_MD", "UL_VLP_RH_MD")]
VLP_covs_MD_l <- melt(pps_VLP_covs_MD, id = c("subj", "group", "sex", "age"), measured = c("LL_VLP_LH_MD", "LL_VLP_RH_MD", "UL_VLP_LH_MD", "UL_VLP_RH_MD"))
names(VLP_covs_MD_l)[6] <- "MD"
VLP_covs_MD_l$hemi <- gl(2, 80, labels = c("L-hemi", "R-hemi"))
VLP_covs_MD_l$limb <- gl(2, 160, labels = c("LOWER", "UPPER"))

VLP_covs_MD_l <- VLP_covs_MD_l[-which(is.na(VLP_covs_MD_l$MD)),]

#model
mod.VLP_grp_limb_covs_MD <- lme(MD ~ -1+group*limb + hemi + sex + age,  random = ~1| subj, data = VLP_covs_MD_l, method = "ML")
summary(mod.VLP_grp_limb_covs_MD)

Linear mixed-effects model fit by maximum likelihood
  Data: VLP_covs_MD_l 
        AIC       BIC   logLik
  -5410.211 -5376.381 2714.106

Random effects:
 Formula: ~1 | subj
         (Intercept)     Residual
StdDev: 2.946014e-05 4.005948e-05

Fixed effects:  MD ~ -1 + group * limb + hemi + sex + age 
                         Value    Std.Error  DF  t-value p-value
group0            0.0007272659 1.782539e-05  76 40.79943  0.0000
group1            0.0007349590 1.782017e-05  76 41.24311  0.0000
limbUPPER        -0.0000094250 6.405072e-06 235 -1.47149  0.1425
hemiR-hemi        0.0000184275 4.552042e-06 235  4.04819  0.0001
sex              -0.0000023262 8.707804e-06  76 -0.26714  0.7901
age               0.0000000493 2.860260e-07  76  0.17244  0.8636
group1:limbUPPER -0.0000098408 9.120401e-06 235 -1.07899  0.2817
 Correlation: 
                 group0 group1 lUPPER hmR-hm sex    age   
group1            0.864                                   
limbUPPER        -0.180  0.000                            
hemiR-hemi       -0.127 -0.125  0.000                     
sex              -0.447 -0.447  0.000 -0.003              
age              -0.862 -0.860  0.000  0.000  0.145       
group1:limbUPPER  0.122 -0.137 -0.702 -0.004  0.000  0.006

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-2.1117077 -0.5613807 -0.1581644  0.3595947  3.9563898 

Number of Observations: 317
Number of Groups: 80