Devoir de stats avancé Hamilton

Reste à faire: graph des distributions, NA clust ?, revoir le modèle mixte. HAMILTON CEST UNE HETEROEVALUATION

1 Importation des tables

library(readxl)
scl <- read_excel("D:/r/m2/hamilton/outils autoeval.xls")
#Nous avons un message R qui nous indique qu'il y a des "ND" dans le data frame (pour non déterminé sans doute), ils sont automatiquement transfomés en NA 
groupe <- read_excel("D:/r/m2/hamilton/outils groupe.xls")
hdrs <- read_excel("D:/r/m2/hamilton/outils hdrs.xls")

2 Data management

On vérifie qu’il n’y ait pas de doublons dans le fichier groupe

table(duplicated(groupe$NUMERO))

## 
## FALSE 
##   146

Pour les fichiers hdrset scl, il est attendu d’avoir des doublons car ce sont des données répétées. En revanche, nous pouvons constater qu’il manque 3 sujets (ou NUMEROS) dans la table scl. Les observations correspondantes seront donc codées en NA lors de la fusion du dataframe complet.

Doublons dans le fichier hdrs:

table(duplicated(hdrs$NUMERO))

## 
## FALSE  TRUE 
##   146   907

Doublons dans le fichier scl:

table(duplicated(scl$NUMERO))

## 
## FALSE  TRUE 
##   143   891

On crée un dataframe dsa, pour “devoir stat avancé”, en fusionnant les tables hdrs et groupe

dsa1<-merge(hdrs,groupe,by = "NUMERO",all.x = T,all.y = T,sort = F)

dsa2<-merge(dsa1,scl,by = c("NUMERO","VISIT"),all.x = T,all.y = T,sort = F)
# scl = 1034 observations vs dsa1 = 1053 obs, on aura donc des NA pour les variables de hdrs pour 19 obs

Par commodité nous changeons l’ordre des colonnes, avec la colonne groupe (n°21) en 2ème colonne

dsa2<-dsa2[,c(1,21,2:20,22:111)]

Repérage d’éventuelles données aberrantes:

Des valeurs >4 sont aberrantes pour les 2 echelles, dont les modalités de réponses vont de 0 à 4. L’échelle HDRS ne comporte pas de données aberrantes sur notre base de donnée, alors que scl90 en comporte de nombreuses. Les numéros de question concernés de l’échelle scl90 sont les suivant:

"%notin%"<-Negate("%in%")

objet1<-list()
for (i in 4:ncol(dsa2)){ 
  if (sum(na.omit(dsa2[,i])%notin%0:4)!=0){
  objet1[i]<-i
  }
}

# Pas de valeurs aberrantes pour l'échelle HDRS
unlist(objet1[4:21]) 

## NULL

#Numéros de question de l'échelle scl90 qui comportent des valeurs aberrantes:
unlist(objet1[-c(1:21)])-21

##  [1]  1  4  5  6  7  8  9 10 12 14 15 16 17 19 20 21 23 24 25 26 28 29 30 32 34
## [26] 37 38 40 41 43 44 46 49 50 51 53 54 55 57 58 63 64 66 68 70 71 73 74 75 76
## [51] 77 79 81 82 83 84 85 86 87 88

Nous pouvons également nous demander si certaines observations ne regroupent pas plusieurs données aberrantes, pour scl90. Voici le nombre de questions dont le score est >4, pour chaque observations:

erreur<-list()
for (i in 22:111){
  erreur[[i]]<-which(dsa2[,i]>4)
  rm(i)
}

sort(table(unlist(erreur,use.names = F)),decreasing = T)

## 
##  364  361  362  134  358  359  360  363    2   70  291  373  379  579  580  582 
##    9    6    4    3    3    3    3    3    2    2    2    2    2    2    2    2 
##  583  841  873  967 1008    3    4   17   21   28   59   66   67  119  121  122 
##    2    2    2    2    2    1    1    1    1    1    1    1    1    1    1    1 
##  126  127  137  239  245  267  290  292  306  332  351  375  376  387  405  423 
##    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
##  463  465  493  494  495  496  498  533  534  548  551  553  555  574  581  584 
##    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
##  586  605  637  642  683  691  697  702  709  718  732  779  802  807  839  862 
##    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
##  875  876  879  881  887  888  889  943  989  993 
##    1    1    1    1    1    1    1    1    1    1

On note par exemple que l’observation 364 comporte des données aberrantes pour 9 items de l’échelle scl90

Nous remplacons ces valeurs aberrantes par des NA.

for (i in 22:111){
  dsa2[dsa2[,i]>4&!is.na(dsa2[,i]),i] <-NA
}

#Vérification en relancant la boucle précédente:
erreur<-list()
for (i in 22:111){
  erreur[[i]]<-which(dsa2[,i]>4)
  rm(i)
}
sort(table(unlist(erreur,use.names = F)),decreasing = T)

## integer(0)

# Remplacement OK

Cas particulier de la question HAMD16A et B de l’échelle HDRS: il fallait remplir l’une ou l’autre. Nous vérifions qu’il n’existe pas d’observations pour lesquelles les 2 sont renseignées:

table(!is.na(dsa2$HAMD16A)&!is.na(dsa2$HAMD16B))

## 
## FALSE 
##  1053

Nous vérifions qu’il n’y ait pas de doublons dans la colonne VISIT pour le même NUMERO

table(duplicated(dsa2[,c("NUMERO","VISIT")]))

## 
## FALSE 
##  1053

On fusionne 16A et 16B, puis on supprime les 2 anciennes variables séparées. On se rend compte qu’il ya une obervation qui ne comporte que des NA : NUMERO 128, grpe 1,J7. Nous ne la supprimons pas.

dsa2$HAMD16<- ifelse(test = is.na(dsa2$HAMD16A),yes = dsa2$HAMD16B,no = dsa2$HAMD16A)

#Vérification
## table(dsa2$HAMD16,exclude=NULL)
## which(is.na(dsa2$HAMD16))

#on supprime les 2 anciennes variables
dsa2_1<-dsa2[,names(dsa2)!="HAMD16A"&names(dsa2)!="HAMD16B"]
# On met ensuite la colonne à la bonne place.
dsa_long<-dsa2_1[,c(1:18,110,19:109)]

3 Présentation des données

Par commodité, et par anticipation de la question suivante, nous utiliserons la base de donnée sous format “large” ou “wide”, c’est à dire avec une ligne par patient:

dsa_wide<-reshape(dsa_long, idvar=c('NUMERO',"GROUPE"), timevar='VISIT', direction='wide')

3.1 Distributions des items de l’échelle HDRS à J0

library(viridis)
library(ggplot2)

p<-list()

for (i in 3:19){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
  scale_fill_viridis(discrete = T,option = "viridis")
    
  )
  p[[i]] <- ggplotGrob(g)
} 

library(gridExtra)

#Necessité de supprimer les éléments NULL de la liste
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))

3.2 Distributions des items de l’échelle HDRS à J56

library(viridis)
library(ggplot2)

p<-list()

for (i in 752:769){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange")

  )
  p[[i]] <- ggplotGrob(g)
} 

library(gridExtra)
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))

3.3 Distributions des items de l’échelle scl90 à J0

#CHUNK A SUPPRIMER
essayer de faire une boucle pour faire tous les 9

seq(0,90,by=9)
 [1]  0  9 18 27 36 45 54 63 72 81 90
vecteurtest<-seq(0,90,by=9)

for (j in vecteurtest){
 for (i in (20:28)+j){
 print(i)
   
 }
  cat("\n")
}
# ça marche, appliquer ça à la boucle en dessous, pour pas à avoir à refaire 10 fois la manip (on a 9 graph par loop, et on a 90 item, donc 10 loop)

#CHUNK A SUPPRIMER
# Ca marche pas, je comrpends pas pourquoi, sachant que ça marche parfaitement pour l'exemple au dessus

library(viridis)
library(ggplot2)

vecteurtest<-seq(0,90,by=9)
 

for (j in vecteurtest){
 p<-list()

 for (i in (20:28)+j){

  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange")
  )
  
  p[[i]] <- ggplotGrob(g)
 
}
library(gridExtra)

marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))
 
}

library(viridis)
library(ggplot2)

 p<-list()
 for (i in 20:108){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange"))
  p[[i]] <- ggplotGrob(g)
}
library(gridExtra)
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))

3.4 Distributions des items de l’échelle scl90 à J56

library(viridis)
library(ggplot2)

 p<-list()
 for (i in 769:858){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange"))
  p[[i]] <- ggplotGrob(g)
}
library(gridExtra)
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))

4 Question 1: Validité de l’echelle de score Hamilton

Nous retrouvons notre base de donnée sous format “large” ou “wide”, créee à la question précédente.

4.0.1 Que vaut la mesure ?

Pour pouvoir faire des sommes globales de score de l’echelle Hamilton, il faut que celle-ci soit unidimensionelle.

Pour évaluer la structure dimensionnelle de l’echelle, on réalise un diagramme des valeurs propres, ou “scree plot”, à partir de la librairie psy

library(psy)

class.source = “fold-show”

D’abord à J0:

scree.plot(dsa_wide[,3:19], simu = 20) # Hamilton J0

Puis à J56:

scree.plot(dsa_wide[,752:768], simu = 20) # Hamilton J56

Nous pouvons remarquer qu’à J0, l’échelle n’est pas du tout unidimensionnelle. On peut compter au moins 3 dimensions, au dessus de l’analyse parallèle.

En revanche à J56, on peut considérer qu’elle est bien unidimensionnelle.

Le score de Hamilton étant une mesure subjective, la fiabilité de l’échelle s’évalue à l’aide du coefficient de corrélation de Cronbach, que l’on peut utiliser en théorie seulement si l’échelle est bien unidimensionnelle.

Cronbach à JO:

cronbach(dsa_wide[,3:19])

## $sample.size
## [1] 146
## 
## $number.of.items
## [1] 17
## 
## $alpha
## [1] 0.45594

Cronbach à J56:

cronbach(dsa_wide[,752:768]) 

## $sample.size
## [1] 120
## 
## $number.of.items
## [1] 17
## 
## $alpha
## [1] 0.819954

4.0.2 Que mesure l’instrument ?

Dans cette section il s’agit d’évaluer la validité concourante et divergente, avec l’echelle SCL90 comme outil de comparaison.

Nous commencons par créer un vecteur du score HDRS contenant le score total pour chaque sujet, d’abord à J0, puis à J56. Puis nous faisons de même avec chaque sous score de l’échelle scl90

4.0.2.1 HDRS et scl90 à J0

Critiquer le fait qu’on puisse faire une somme total pour hdrs, en regard de l’allure du screeplot

Pour scl90, nous remettons ici les items correspondants à chaque sous scores:

depression 5,14,15,20,22,26,29,30,31,32,54,71,79
Somatisation 1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40
Symptômes obsessionnels 9, 10, 28, 38, 3,45, 46, 51, 55, 65
Sensitivité interpersonnelle ou vulnérabilité 6, 21, 34, 36, 37,41, 61, 69, 73
Anxiété 2, 17, 23, 33, 39,57, 72, 78, 80, 86
Hostilité 11, 24, 63, 67,74, 81
Phobies 13, 25, 47, 70,75, 82, 50
Traits paranoïaques 8, 18, 43, 68, 76,83
Traits psychotiques 7, 16, 35, 62, 77, 84, 85, 87, 90, 88
Symptômes divers 19, 44, 59, 60, 64, 66, 89

Pour ensuite calculer les sous score HDRS et scl90

# HDRS J0:
score_tot_hdrs_j0<-c()

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,3:19]))==0){
    score_tot_hdrs_j0[i]<-sum(dsa_wide[i,3:19])
  } else {
    score_tot_hdrs_j0[i]<-NA
  }
}



score_tot_scl_dep_j0<-c()
score_tot_scl_soma_j0<-c()
score_tot_scl_obses_j0<-c()
score_tot_scl_sensit_j0<-c()
score_tot_scl_anx_j0<-c()
score_tot_scl_host_j0<-c()
score_tot_scl_phob_j0<-c()
score_tot_scl_paran_j0<-c()
score_tot_scl_psycho_j0<-c()
score_tot_scl_divers_j0<-c()


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+19]))==0){
    score_tot_scl_dep_j0[i]<-sum(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+19])
  }else{
    score_tot_scl_dep_j0[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+19]))==0){
    score_tot_scl_soma_j0[i]<-sum(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+19])
  }else{
    score_tot_scl_soma_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+19]))==0){
    score_tot_scl_obses_j0[i]<-sum(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+19])
  }else{
    score_tot_scl_obses_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+19]))==0){
    score_tot_scl_sensit_j0[i]<-sum(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+19])
  }else{
    score_tot_scl_sensit_j0[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19]))==0){
    score_tot_scl_anx_j0[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19])
  }else{
    score_tot_scl_anx_j0[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19]))==0){
    score_tot_scl_anx_j0[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19])
  }else{
    score_tot_scl_anx_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(11, 24, 63, 67,74, 81)+19]))==0){
    score_tot_scl_host_j0[i]<-sum(dsa_wide[i,c(11, 24, 63, 67,74, 81)+19])
  }else{
    score_tot_scl_host_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+19]))==0){
    score_tot_scl_phob_j0[i]<-sum(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+19])
  }else{
    score_tot_scl_phob_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(8, 18, 43, 68, 76,83)+19]))==0){
    score_tot_scl_paran_j0[i]<-sum(dsa_wide[i,c(8, 18, 43, 68, 76,83)+19])
  }else{
    score_tot_scl_paran_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+19]))==0){
    score_tot_scl_psycho_j0[i]<-sum(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+19])
  }else{
    score_tot_scl_psycho_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+19]))==0){
    score_tot_scl_divers_j0[i]<-sum(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+19])
  }else{
    score_tot_scl_divers_j0[i]<-NA
  }
}

#Nous créons un data frame qui rassemble l'ensemble de ces vecteurs : "df_score_tot_scl_j0"

df_score_tot_scl_j0<-data.frame(
  score_tot_scl_dep_j0,
  score_tot_scl_soma_j0,
  score_tot_scl_obses_j0,
  score_tot_scl_sensit_j0,
  score_tot_scl_anx_j0,
  score_tot_scl_host_j0,
  score_tot_scl_phob_j0,
  score_tot_scl_paran_j0,
  score_tot_scl_psycho_j0,
  score_tot_scl_divers_j0
  )

Nous pouvons ensuite calculer la correlation du score HDRS à J0 avec l’ensemble des sous score scl90 à J0. Nous obtenons le tableau de corrélations suivant:

table_corr_j0<-data.frame()
table_col_names_j0<-c()

for (i in 1:ncol(df_score_tot_scl_j0)){
  table_corr_j0[1,i]<-round(cor(score_tot_hdrs_j0,df_score_tot_scl_j0[,i],use = "complete.obs"),digits = 3)
  table_col_names_j0[i]<-names(df_score_tot_scl_j0[i])
}

rm(i)

names(table_corr_j0)<-table_col_names_j0

library(knitr)
kable(t(table_corr_j0),caption = "Corrélations entre le score HDRS et tous les sous score de l'échelle scl90, à J0")

Corrélations entre le score HDRS et tous les sous score de l’échelle scl90, à J0

	1
score_tot_scl_dep_j0	0.278
score_tot_scl_soma_j0	0.367
score_tot_scl_obses_j0	0.234
score_tot_scl_sensit_j0	-0.105
score_tot_scl_anx_j0	0.135
score_tot_scl_host_j0	0.060
score_tot_scl_phob_j0	0.075
score_tot_scl_paran_j0	-0.112
score_tot_scl_psycho_j0	-0.029
score_tot_scl_divers_j0	0.432

4.0.2.2 HDRS et scl90 à J56

En procédant de la même manière pour les scores à J56, nous obtenons le tableau de corrélations suivant:

# score hdrs J56:
score_tot_hdrs_j56<-c()

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,752:768]))==0){
    score_tot_hdrs_j56[i]<-sum(dsa_wide[i,752:768])
  } else {
    score_tot_hdrs_j56[i]<-NA
  }
}


# Sous scores scl90 à J56:

score_tot_scl_dep_j56<-c()
score_tot_scl_soma_j56<-c()
score_tot_scl_obses_j56<-c()
score_tot_scl_sensit_j56<-c()
score_tot_scl_anx_j56<-c()
score_tot_scl_host_j56<-c()
score_tot_scl_phob_j56<-c()
score_tot_scl_paran_j56<-c()
score_tot_scl_psycho_j56<-c()
score_tot_scl_divers_j56<-c()

# pour adapter au dataframe scl j56 on fait +768 
for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+768]))==0){
    score_tot_scl_dep_j56[i]<-sum(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+768])
  }else{
    score_tot_scl_dep_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+768]))==0){
    score_tot_scl_soma_j56[i]<-sum(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+768])
  }else{
    score_tot_scl_soma_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+768]))==0){
    score_tot_scl_obses_j56[i]<-sum(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+768])
  }else{
    score_tot_scl_obses_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+768]))==0){
    score_tot_scl_sensit_j56[i]<-sum(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+768])
  }else{
    score_tot_scl_sensit_j56[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768]))==0){
    score_tot_scl_anx_j56[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768])
  }else{
    score_tot_scl_anx_j56[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768]))==0){
    score_tot_scl_anx_j56[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768])
  }else{
    score_tot_scl_anx_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(11, 24, 63, 67,74, 81)+768]))==0){
    score_tot_scl_host_j56[i]<-sum(dsa_wide[i,c(11, 24, 63, 67,74, 81)+768])
  }else{
    score_tot_scl_host_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+768]))==0){
    score_tot_scl_phob_j56[i]<-sum(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+768])
  }else{
    score_tot_scl_phob_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(8, 18, 43, 68, 76,83)+768]))==0){
    score_tot_scl_paran_j56[i]<-sum(dsa_wide[i,c(8, 18, 43, 68, 76,83)+768])
  }else{
    score_tot_scl_paran_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+768]))==0){
    score_tot_scl_psycho_j56[i]<-sum(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+768])
  }else{
    score_tot_scl_psycho_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+768]))==0){
    score_tot_scl_divers_j56[i]<-sum(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+768])
  }else{
    score_tot_scl_divers_j56[i]<-NA
  }
}

df_score_tot_scl_j56<-data.frame(
  score_tot_scl_dep_j56,
  score_tot_scl_soma_j56,
  score_tot_scl_obses_j56,
  score_tot_scl_sensit_j56,
  score_tot_scl_anx_j56,
  score_tot_scl_host_j56,
  score_tot_scl_phob_j56,
  score_tot_scl_paran_j56,
  score_tot_scl_psycho_j56,
  score_tot_scl_divers_j56
)

table_corr_j56<-data.frame()
table_col_names_j56<-c()

for (i in 1:ncol(df_score_tot_scl_j56)){
  table_corr_j56[1,i]<-round(cor(score_tot_hdrs_j56,df_score_tot_scl_j56[,i],use = "complete.obs"),digits = 3)
  table_col_names_j56[i]<-names(df_score_tot_scl_j56[i])
}

rm(i)

names(table_corr_j56)<-table_col_names_j56


library(knitr)
kable(t(table_corr_j56),caption = "Corrélations entre le score HDRS et tous les sous score de l'échelle scl90, à J56")

Corrélations entre le score HDRS et tous les sous score de l’échelle scl90, à J56

	1
score_tot_scl_dep_j56	0.654
score_tot_scl_soma_j56	0.617
score_tot_scl_obses_j56	0.646
score_tot_scl_sensit_j56	0.499
score_tot_scl_anx_j56	0.494
score_tot_scl_host_j56	0.270
score_tot_scl_phob_j56	0.446
score_tot_scl_paran_j56	0.393
score_tot_scl_psycho_j56	0.531
score_tot_scl_divers_j56	0.633

5 Graphiques

6 Question 2

6.1 LOCF

Il y a des numéros qui sont dispersés et qui ne sont pas à la suite, on va donc les ordonner. On fait de même pour l’ordre des VISIT pour chaque numéro.

#il y a des numéros dispersés
for( i in unique(dsa_long[,1])){
  print(which(dsa_long[,1]==i))
}

#on doit d'abord créer un autre dataframe que l'on va ensuite pouvoir ordonner correctement

#Pour pouvoir ordonner, le plus simple ici est de supprimer le premier character "J"

dsa_long_visit_numeric <-data.frame(dsa_long)
dsa_long_visit_numeric$VISIT<-as.numeric(substring(dsa_long_visit_numeric$VISIT,first = 2))


#on ordonne le tout:
dsa_long_sorted_by_numero_and_visit<-dsa_long_visit_numeric[order(dsa_long_visit_numeric$NUMERO,dsa_long_visit_numeric$VISIT),]

On vérifie que les opérations se soient correctement déroulées

# on vérifie:
unique(dsa_long_sorted_by_numero_and_visit$NUMERO)
unique(dsa_long_sorted_by_numero_and_visit$VISIT)

# on vérifie aussi qu'on n'a plus de numéros dispérsés parmi les lignes:
for( i in unique(dsa_long_sorted_by_numero_and_visit[,1])){
  print(which(dsa_long_sorted_by_numero_and_visit[,1]==i))
}

# et qu'on a bien toutes nos observ:
dim(dsa_long_sorted_by_numero_and_visit)

On réinitialise le numéro des lignes pour plus de clarté.

dsa_long_sorted_by_numero_and_visit_rebootrowname<- dsa_long_sorted_by_numero_and_visit
rownames(dsa_long_sorted_by_numero_and_visit_rebootrowname)<-NULL

Dans les données disponibles, les visites non effectuées ne sont pas renseignées, nous n’avons donc pas de “NA” pour ces visites, excépté pour une observation. Nous pouvons donc prendre simplement la dernière visite renseignée pour chaque numéros.

#d'abord on renome le df pour plus de lisibilité:
dsa_long_final<-dsa_long_sorted_by_numero_and_visit_rebootrowname

#Nous allons avoir recours à une liste, nous décidons de changer le numéro "0" pour "1"
dsa_long_final2<-dsa_long_final
dsa_long_final2[1,1]<-1


# Nous indexons les scores pour chaque numéros et chaques visites 

index_locf<-list()
for (i in unique(dsa_long_final2$NUMERO)){
  index_locf[[i]]<-unique(dsa_long_final2[dsa_long_final2$NUMERO==i,3])

}

#pour pouvoir ensuite calculer les sommes totales 




tot_locf <- vector(mode="list", length=length(index_locf))
for (i in 1:length(index_locf)){
for (j in 1:length(index_locf[[i]])){
 tot_locf[[i]][j] <-sum(dsa_long_final2[dsa_long_final2$NUMERO==i&dsa_long_final2$VISIT==index_locf[[i]][j],4:20])
}
}


#et prendre le score correspondant à la dernière visite. (On utilisera na.locf pour le NA à la dernière visite du numéro 128)

"%notin%"<-Negate("%in%")
library(zoo)

tot_locf_done_j56<-c()
for (i in 1:length(tot_locf)){
  if(i %notin% which(sapply(index_locf,is.null))){
  tot_locf_done_j56[i]<-tail(na.locf(unlist(tot_locf[i])),1)
  }
}
tot_locf_done_hdrs_j56<-na.omit(tot_locf_done_j56)

On obtient un vecteur dans l’ordre des NUMERO triés par ordre croissant

Nous allons donc appliquer cet ordre au fichier large dénormalisé (dsa_wide), et également changer le numéro 0 en numéro 1, avant de fusionner le vecteur pour créer une nouvelle colonne.

dsa_wide_sorted_by_numero<-dsa_wide[order(dsa_wide$NUMERO),]

# on remplace le numero 0 par 1:
dsa_wide_sorted_by_numero_0removed<-dsa_wide_sorted_by_numero
dsa_wide_sorted_by_numero_0removed[1,1]<-1

#on réinitialise là aussi les numéros de lignes:
dsa_wide_sorted_by_numero_0removed_rebootrowname<-dsa_wide_sorted_by_numero_0removed
rownames(dsa_wide_sorted_by_numero_0removed_rebootrowname)<-NULL

score_tot_hdrs_j0_sorted<-c()
for (i in unique(dsa_long_final2$NUMERO)){
score_tot_hdrs_j0_sorted[i]<-sum(dsa_long_final2[dsa_long_final2$NUMERO==i&dsa_long_final2$VISIT==0,4:20])
}
score_tot_hdrs_j0_sorted_naomit<-na.omit(score_tot_hdrs_j0_sorted)

Nous pouvons enfin fusionner au fichier long les vecteurs des scores HDRS à J0 et J56

dsa_wide_merged<-cbind(dsa_wide_sorted_by_numero_0removed_rebootrowname,score_tot_hdrs_j0_sorted_naomit,tot_locf_done_hdrs_j56)

Avec la fonction ‘unique()’ qui permet d’afficher des résultats par ordre d’occurence, nous pouvons vérifier que nous avons bien le même ordre de numéros dans le fichier long et large, que les numéros “match”

check<-c()
for (i in 1:length(dsa_wide_sorted_by_numero_0removed_rebootrowname$NUMERO)){
  check[i]<-unique(dsa_wide_sorted_by_numero_0removed_rebootrowname$NUMERO)[i]==unique(dsa_long_final2$NUMERO)[i]
}

kable(table(check), caption="Match")

Match

check	Freq
TRUE	146

Avec la méthode LOCF, nous considérons qu’il n’y a pas de perdus de vu, on peut donc procéder à une simple comparaison des moyennes des score HDRS entre les 2 groupes, à J0 et J56

ATTENTION l’histogramme de score_tot_hdrs_j56 est très asymétrique, on considère la robustesse lorsque n>30 ou bien on fait un bootstrap ? Ca à un sens de faire un bootstrap ici ?

Test t de student à J0:

t.test(dsa_wide_merged[dsa_wide_merged[,2]==0,859],dsa_wide_merged[dsa_wide_merged[,2]==1,859],var.equal = T)

## 
##  Two Sample t-test
## 
## data:  dsa_wide_merged[dsa_wide_merged[, 2] == 0, 859] and dsa_wide_merged[dsa_wide_merged[, 2] == 1, 859]
## t = 0.37947, df = 144, p-value = 0.7049
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.113645  1.642847
## sample estimates:
## mean of x mean of y 
##  28.05333  27.78873

Test t de student à J56:

t.test(dsa_wide_merged[dsa_wide_merged[,2]==0,860],dsa_wide_merged[dsa_wide_merged[,2]==1,860],var.equal = T)

## 
##  Two Sample t-test
## 
## data:  dsa_wide_merged[dsa_wide_merged[, 2] == 0, 860] and dsa_wide_merged[dsa_wide_merged[, 2] == 1, 860]
## t = -3.7275, df = 144, p-value = 0.0002772
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.447292 -2.286042
## sample estimates:
## mean of x mean of y 
##  8.133333 13.000000

6.2 Modèle mixte

Dans un modèle de régression pour expliquer le score HDRS, nous devons ici introduire la varibale NUMERO comme variable explicative car ce sont des données répétées et donc corrélées. Nous pouvons coder cette variable en qualitative nominale (avec donc k-1 variable binaire dans le modèle) ou bien la coder en effet aléatoire, ce qui est conseillé lorsqu’on a plus de 15/20 modalités (ce qui est le cas ici avec 146 NUMEROS)

Il nous faut utiliser le fichier long car nous avons besoin de la variable VISIT. Nous aurons également besoin du score total HDRS pour chaque NUMEROS, à chaque VISIT, que nous rajouterons dans une nouvelle colonne (à la fin pour pas déranger l’ordre des colonnes). Nous avons déjà J0 et J56 calculé précédemment

tot_par_visit<-c()
for (i in 1:nrow(dsa_long_final2)){
  tot_par_visit[i]<-sum(dsa_long_final2[i,4:20],na.rm = F)
}

## length(tot_par_visit) # on a bien 1053 valeur qui correspond au nombre d'observations, dans l'ordre des lignes, en ayant gardé le NA

# on ajoute ce vecteur en tant que dernière colonne du Df:
dsa_long_final_tot_par_visit<-cbind(dsa_long_final2,tot_par_visit)

Pour utiliser un modèle mixte, nous aurons besoin de transformer VISIT en facteur, puis en numérique.

dsa_long_final_tot_par_visit_factor<-data.frame(dsa_long_final_tot_par_visit)

dsa_long_final_tot_par_visit_factor$VISIT<-as.numeric(as.factor(dsa_long_final_tot_par_visit_factor$VISIT))

#on se rassure: il n'ya pas eu de changement d'ordre dans cette manip, exemple: as.numeric(as.factor(c(1,98,32))) = [1] 1 3 2

On utilisera le package lmerTest, qui nous donne les p values des tests

library(lmerTest)

summary(lmer(formula = tot_par_visit~GROUPE+VISIT+(VISIT|NUMERO), data = dsa_long_final_tot_par_visit_factor))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: tot_par_visit ~ GROUPE + VISIT + (VISIT | NUMERO)
##    Data: dsa_long_final_tot_par_visit_factor
## 
## REML criterion at convergence: 6283.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0382 -0.5588  0.0112  0.5359  4.3903 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  NUMERO   (Intercept) 25.6312  5.0627        
##           VISIT        0.6612  0.8132   -0.48
##  Residual             14.6416  3.8264        
## Number of obs: 1052, groups:  NUMERO, 146
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  28.97343    0.62486 165.85293  46.368   <2e-16 ***
## GROUPE        2.21152    0.78938 141.81953   2.802   0.0058 ** 
## VISIT        -2.94911    0.08833 134.10575 -33.387   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) GROUPE
## GROUPE -0.616       
## VISIT  -0.478  0.009

On obtient une p value significative à 0.0058 pour l’effet fixe GROUPE avec un coeff à 2.21152, ce qui signifie que lorsqu’on passe du grpe 0 à 1 le score HDRS augmente de 2.21152.

On peut remarquer que cette différence est comprise dans l’intervalle de confiance du calcul du t test, réalisé avec locf (-7.447292 -2.286042)

L’utilisation de ce modèle repose sur une hypothèse de linéarité de l’évolution des scores, du fait des perdus de vu. Voyons ce qu’il en est. Nous présentons ci-dessous l’évolution des scores pour un échantillon de 30 patients tirés au sort parmi les NUMERO du dataframe, afin d’apprécier la linéarité de cette évolution:

library(ggplot2)
library(viridis)

ggplot(data=dsa_long_final_tot_par_visit_factor[dsa_long_final_tot_par_visit_factor$NUMERO %in%
                                                  unique(dsa_long_final_tot_par_visit_factor$NUMERO)[sample(x = unique(dsa_long_final_tot_par_visit_factor$NUMERO),size = 30,replace = F)],],
       mapping = aes(x = VISIT,y =tot_par_visit, group=as.factor(NUMERO),color=as.factor(NUMERO)) )+
  geom_line()+
  scale_colour_viridis(discrete = T)+
  theme(legend.position = "none")

Voici le même graphique avec représentation de tous les sujets par groupes, et la courbe des scores moyens en fonction du temps de visite

#Il nous faut créer un nouveau dataframne, pour ajouter un geom line spécifique

mean_group_visit<-aggregate(tot_par_visit~GROUPE+VISIT,dsa_long_final_tot_par_visit_factor,mean)


ggplot(data =dsa_long_final_tot_par_visit_factor[dsa_long_final_tot_par_visit_factor$NUMERO%in%
                                                   unique(dsa_long_final_tot_par_visit_factor$NUMERO),],
       mapping = aes(x = VISIT,y = tot_par_visit, color=as.factor(GROUPE)) )+
  geom_line(aes(group=as.factor(NUMERO)),alpha=0.4)+
  geom_line(data = mean_group_visit,size=3)+
  scale_colour_viridis(discrete = T)+
  theme(legend.position = "none")

7 Question 3

Nous choisirons de coder une diminution de 50% du score HDRS en 1, et 0 sinon. Nous stockerons les résultats dans une nouvelle colonne.

cbinaire<-vector(mode="list", length=191)

for ( i in dsa_long_final_tot_par_visit_factor$NUMERO){
  for(j in dsa_long_final_tot_par_visit_factor[
    dsa_long_final_tot_par_visit_factor$NUMERO==i,3]){
    
      
   if(is.na(dsa_long_final_tot_par_visit_factor[
        dsa_long_final_tot_par_visit_factor$NUMERO==i&
          dsa_long_final_tot_par_visit_factor$VISIT==j,111])){
        cbinaire[[i]][j]<-NA
   }
    else if 
      (dsa_long_final_tot_par_visit_factor[
        dsa_long_final_tot_par_visit_factor$NUMERO==i&
        dsa_long_final_tot_par_visit_factor$VISIT==j,111] >
       dsa_long_final_tot_par_visit_factor[
         dsa_long_final_tot_par_visit_factor$NUMERO==i&
         dsa_long_final_tot_par_visit_factor$VISIT==1,111]/2){
        cbinaire[[i]][j]<-0
    }
     else{
       cbinaire[[i]][j]<-1
     }
  }
}


# Pour une raison indéterminée, 2 NA se sont glissés en plus dans les données, ce qui ramène les observations à 1055 au lieu de 1053. On les supprime:
## which(sapply(X = cbinaire,FUN = function(X) NA %in% X))


cbinaire_f<-cbinaire
cbinaire_f[[92]]<-cbinaire_f[[92]][-2]
cbinaire_f[[170]]<-cbinaire_f[[170]][-3]


# Les résultats sont stockés en dernière colonne du data frame:

dsa_long_final_tot_par_visit_factor_2<-cbind(dsa_long_final_tot_par_visit_factor,unlist(cbinaire_f))

Nous pouvons visulaiser l’évolution des deux groupes à l’aide de la méthode de Kaplan Meier, et réaliser un test du log-rank

library(survival)

#Kaplan Meier
plot(survfit(Surv(VISIT,unlist(cbinaire_f))~GROUPE,data=dsa_long_final_tot_par_visit_factor_2),col=c("red","blue"))

#test du log-rank
survdiff(Surv(VISIT,unlist(cbinaire_f))~GROUPE,data=dsa_long_final_tot_par_visit_factor_2)

## Call:
## survdiff(formula = Surv(VISIT, unlist(cbinaire_f)) ~ GROUPE, 
##     data = dsa_long_final_tot_par_visit_factor_2)
## 
## n=1052, 1 observation effacée parce que manquante.
## 
##            N Observed Expected (O-E)^2/E (O-E)^2/V
## GROUPE=0 554      273      229      8.33      29.2
## GROUPE=1 498      157      201      9.52      29.2
## 
##  Chisq= 29.2  on 1 degrees of freedom, p= 6e-08

8 Appendice

knitr::opts_chunk$set(warning=FALSE, message=FALSE)
library(readxl)
scl <- read_excel("D:/r/m2/hamilton/outils autoeval.xls")
#Nous avons un message R qui nous indique qu'il y a des "ND" dans le data frame (pour non déterminé sans doute), ils sont automatiquement transfomés en NA 
groupe <- read_excel("D:/r/m2/hamilton/outils groupe.xls")
hdrs <- read_excel("D:/r/m2/hamilton/outils hdrs.xls")
table(duplicated(groupe$NUMERO))
table(duplicated(hdrs$NUMERO))

table(duplicated(scl$NUMERO))
dsa1<-merge(hdrs,groupe,by = "NUMERO",all.x = T,all.y = T,sort = F)
dsa2<-merge(dsa1,scl,by = c("NUMERO","VISIT"),all.x = T,all.y = T,sort = F)
# scl = 1034 observations vs dsa1 = 1053 obs, on aura donc des NA pour les variables de hdrs pour 19 obs
dsa2<-dsa2[,c(1,21,2:20,22:111)]
"%notin%"<-Negate("%in%")

objet1<-list()
for (i in 4:ncol(dsa2)){ 
  if (sum(na.omit(dsa2[,i])%notin%0:4)!=0){
  objet1[i]<-i
  }
}

# Pas de valeurs aberrantes pour l'échelle HDRS
unlist(objet1[4:21]) 

#Numéros de question de l'échelle scl90 qui comportent des valeurs aberrantes:
unlist(objet1[-c(1:21)])-21

erreur<-list()
for (i in 22:111){
  erreur[[i]]<-which(dsa2[,i]>4)
  rm(i)
}

sort(table(unlist(erreur,use.names = F)),decreasing = T)

for (i in 22:111){
  dsa2[dsa2[,i]>4&!is.na(dsa2[,i]),i] <-NA
}

#Vérification en relancant la boucle précédente:
erreur<-list()
for (i in 22:111){
  erreur[[i]]<-which(dsa2[,i]>4)
  rm(i)
}
sort(table(unlist(erreur,use.names = F)),decreasing = T)
# Remplacement OK
table(!is.na(dsa2$HAMD16A)&!is.na(dsa2$HAMD16B))
table(duplicated(dsa2[,c("NUMERO","VISIT")]))

dsa2$HAMD16<- ifelse(test = is.na(dsa2$HAMD16A),yes = dsa2$HAMD16B,no = dsa2$HAMD16A)

#Vérification
table(dsa2$HAMD16,exclude=NULL)
which(is.na(dsa2$HAMD16))

#on supprime les 2 anciennes variables 
dsa2_1<-dsa2[,names(dsa2)!="HAMD16A"&names(dsa2)!="HAMD16B"]
# On met ensuite la colonne à la bonne place.
dsa_long<-dsa2_1[,c(1:18,110,19:109)]
dsa_wide<-reshape(dsa_long, idvar=c('NUMERO',"GROUPE"), timevar='VISIT', direction='wide')
library(viridis)
library(ggplot2)

p<-list()

for (i in 3:19){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
  scale_fill_viridis(discrete = T,option = "viridis")
    
  )
  p[[i]] <- ggplotGrob(g)
} 

library(gridExtra)

#Necessité de supprimer les éléments NULL de la liste
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))

library(viridis)
library(ggplot2)

p<-list()

for (i in 752:769){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange")

  )
  p[[i]] <- ggplotGrob(g)
} 

library(gridExtra)
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))
#CHUNK A SUPPRIMER
essayer de faire une boucle pour faire tous les 9

seq(0,90,by=9)
 [1]  0  9 18 27 36 45 54 63 72 81 90
vecteurtest<-seq(0,90,by=9)

for (j in vecteurtest){
 for (i in (20:28)+j){
 print(i)
   
 }
  cat("\n")
}
# ça marche, appliquer ça à la boucle en dessous, pour pas à avoir à refaire 10 fois la manip (on a 9 graph par loop, et on a 90 item, donc 10 loop)
#CHUNK A SUPPRIMER
# Ca marche pas, je comrpends pas pourquoi, sachant que ça marche parfaitement pour l'exemple au dessus

library(viridis)
library(ggplot2)

vecteurtest<-seq(0,90,by=9)
 

for (j in vecteurtest){
 p<-list()

 for (i in (20:28)+j){

  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange")
  )
  
  p[[i]] <- ggplotGrob(g)
 
}
library(gridExtra)

marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))
 
}

library(viridis)
library(ggplot2)

 p<-list()
 for (i in 20:108){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange"))
  p[[i]] <- ggplotGrob(g)
}
library(gridExtra)
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))


library(viridis)
library(ggplot2)

 p<-list()
 for (i in 769:858){
  g<-(ggplot(dsa_wide, aes(x='', (..count..) / sum(..count..), fill=as.factor(dsa_wide[,i]))) + 
geom_bar(position = position_dodge())+
  labs(x=names(dsa_wide[i]))+
  theme(axis.title.y = element_blank())+
   theme(legend.title =element_blank())+
 theme(legend.key.size = unit(0.5, 'cm'))+
    scale_fill_viridis(discrete = T,option = "viridis",na.value="orange"))
  p[[i]] <- ggplotGrob(g)
}
library(gridExtra)
marrangeGrob(grobs=p[lengths(p) != 0], layout_matrix=matrix(c(1:9),nrow = 3,ncol = 3,byrow = T))

library(psy)
scree.plot(dsa_wide[,3:19], simu = 20) # Hamilton J0
scree.plot(dsa_wide[,752:768], simu = 20) # Hamilton J56

cronbach(dsa_wide[,3:19])
cronbach(dsa_wide[,752:768]) 
# HDRS J0:
score_tot_hdrs_j0<-c()

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,3:19]))==0){
    score_tot_hdrs_j0[i]<-sum(dsa_wide[i,3:19])
  } else {
    score_tot_hdrs_j0[i]<-NA
  }
}



score_tot_scl_dep_j0<-c()
score_tot_scl_soma_j0<-c()
score_tot_scl_obses_j0<-c()
score_tot_scl_sensit_j0<-c()
score_tot_scl_anx_j0<-c()
score_tot_scl_host_j0<-c()
score_tot_scl_phob_j0<-c()
score_tot_scl_paran_j0<-c()
score_tot_scl_psycho_j0<-c()
score_tot_scl_divers_j0<-c()


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+19]))==0){
    score_tot_scl_dep_j0[i]<-sum(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+19])
  }else{
    score_tot_scl_dep_j0[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+19]))==0){
    score_tot_scl_soma_j0[i]<-sum(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+19])
  }else{
    score_tot_scl_soma_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+19]))==0){
    score_tot_scl_obses_j0[i]<-sum(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+19])
  }else{
    score_tot_scl_obses_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+19]))==0){
    score_tot_scl_sensit_j0[i]<-sum(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+19])
  }else{
    score_tot_scl_sensit_j0[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19]))==0){
    score_tot_scl_anx_j0[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19])
  }else{
    score_tot_scl_anx_j0[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19]))==0){
    score_tot_scl_anx_j0[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+19])
  }else{
    score_tot_scl_anx_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(11, 24, 63, 67,74, 81)+19]))==0){
    score_tot_scl_host_j0[i]<-sum(dsa_wide[i,c(11, 24, 63, 67,74, 81)+19])
  }else{
    score_tot_scl_host_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+19]))==0){
    score_tot_scl_phob_j0[i]<-sum(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+19])
  }else{
    score_tot_scl_phob_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(8, 18, 43, 68, 76,83)+19]))==0){
    score_tot_scl_paran_j0[i]<-sum(dsa_wide[i,c(8, 18, 43, 68, 76,83)+19])
  }else{
    score_tot_scl_paran_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+19]))==0){
    score_tot_scl_psycho_j0[i]<-sum(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+19])
  }else{
    score_tot_scl_psycho_j0[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+19]))==0){
    score_tot_scl_divers_j0[i]<-sum(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+19])
  }else{
    score_tot_scl_divers_j0[i]<-NA
  }
}

#Nous créons un data frame qui rassemble l'ensemble de ces vecteurs : "df_score_tot_scl_j0"

df_score_tot_scl_j0<-data.frame(
  score_tot_scl_dep_j0,
  score_tot_scl_soma_j0,
  score_tot_scl_obses_j0,
  score_tot_scl_sensit_j0,
  score_tot_scl_anx_j0,
  score_tot_scl_host_j0,
  score_tot_scl_phob_j0,
  score_tot_scl_paran_j0,
  score_tot_scl_psycho_j0,
  score_tot_scl_divers_j0
  )

table_corr_j0<-data.frame()
table_col_names_j0<-c()

for (i in 1:ncol(df_score_tot_scl_j0)){
  table_corr_j0[1,i]<-round(cor(score_tot_hdrs_j0,df_score_tot_scl_j0[,i],use = "complete.obs"),digits = 3)
  table_col_names_j0[i]<-names(df_score_tot_scl_j0[i])
}

rm(i)

names(table_corr_j0)<-table_col_names_j0

library(knitr)
kable(t(table_corr_j0),caption = "Corrélations entre le score HDRS et tous les sous score de l'échelle scl90, à J0")
# score hdrs J56:
score_tot_hdrs_j56<-c()

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,752:768]))==0){
    score_tot_hdrs_j56[i]<-sum(dsa_wide[i,752:768])
  } else {
    score_tot_hdrs_j56[i]<-NA
  }
}


# Sous scores scl90 à J56:

score_tot_scl_dep_j56<-c()
score_tot_scl_soma_j56<-c()
score_tot_scl_obses_j56<-c()
score_tot_scl_sensit_j56<-c()
score_tot_scl_anx_j56<-c()
score_tot_scl_host_j56<-c()
score_tot_scl_phob_j56<-c()
score_tot_scl_paran_j56<-c()
score_tot_scl_psycho_j56<-c()
score_tot_scl_divers_j56<-c()

# pour adapter au dataframe scl j56 on fait +768 
for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+768]))==0){
    score_tot_scl_dep_j56[i]<-sum(dsa_wide[i,c(5,14,15,20,22,26,29,30,31,32,54,71,79)+768])
  }else{
    score_tot_scl_dep_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+768]))==0){
    score_tot_scl_soma_j56[i]<-sum(dsa_wide[i,c(1, 4, 12, 27, 42,48, 49, 52, 53,56, 58, 40)+768])
  }else{
    score_tot_scl_soma_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+768]))==0){
    score_tot_scl_obses_j56[i]<-sum(dsa_wide[i,c(9, 10, 28, 38, 3,45, 46, 51, 55, 65)+768])
  }else{
    score_tot_scl_obses_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+768]))==0){
    score_tot_scl_sensit_j56[i]<-sum(dsa_wide[i,c(6, 21, 34, 36, 37,41, 61, 69, 73)+768])
  }else{
    score_tot_scl_sensit_j56[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768]))==0){
    score_tot_scl_anx_j56[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768])
  }else{
    score_tot_scl_anx_j56[i]<-NA
  }
}


for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768]))==0){
    score_tot_scl_anx_j56[i]<-sum(dsa_wide[i,c(2, 17, 23, 33, 39,57, 72, 78, 80, 86)+768])
  }else{
    score_tot_scl_anx_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(11, 24, 63, 67,74, 81)+768]))==0){
    score_tot_scl_host_j56[i]<-sum(dsa_wide[i,c(11, 24, 63, 67,74, 81)+768])
  }else{
    score_tot_scl_host_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+768]))==0){
    score_tot_scl_phob_j56[i]<-sum(dsa_wide[i,c(13, 25, 47, 70,75, 82, 50)+768])
  }else{
    score_tot_scl_phob_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(8, 18, 43, 68, 76,83)+768]))==0){
    score_tot_scl_paran_j56[i]<-sum(dsa_wide[i,c(8, 18, 43, 68, 76,83)+768])
  }else{
    score_tot_scl_paran_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+768]))==0){
    score_tot_scl_psycho_j56[i]<-sum(dsa_wide[i,c( 7, 16, 35, 62, 77, 84, 85, 87, 90, 88)+768])
  }else{
    score_tot_scl_psycho_j56[i]<-NA
  }
}

for (i in 1:nrow(dsa_wide)){
  if (sum(is.na(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+768]))==0){
    score_tot_scl_divers_j56[i]<-sum(dsa_wide[i,c(19, 44, 59, 60, 64, 66, 89)+768])
  }else{
    score_tot_scl_divers_j56[i]<-NA
  }
}

df_score_tot_scl_j56<-data.frame(
  score_tot_scl_dep_j56,
  score_tot_scl_soma_j56,
  score_tot_scl_obses_j56,
  score_tot_scl_sensit_j56,
  score_tot_scl_anx_j56,
  score_tot_scl_host_j56,
  score_tot_scl_phob_j56,
  score_tot_scl_paran_j56,
  score_tot_scl_psycho_j56,
  score_tot_scl_divers_j56
)

table_corr_j56<-data.frame()
table_col_names_j56<-c()

for (i in 1:ncol(df_score_tot_scl_j56)){
  table_corr_j56[1,i]<-round(cor(score_tot_hdrs_j56,df_score_tot_scl_j56[,i],use = "complete.obs"),digits = 3)
  table_col_names_j56[i]<-names(df_score_tot_scl_j56[i])
}

rm(i)

names(table_corr_j56)<-table_col_names_j56


library(knitr)
kable(t(table_corr_j56),caption = "Corrélations entre le score HDRS et tous les sous score de l'échelle scl90, à J56")
#il y a des numéros dispersés
for( i in unique(dsa_long[,1])){
  print(which(dsa_long[,1]==i))
}
#on doit d'abord créer un autre dataframe que l'on va ensuite pouvoir ordonner correctement

#Pour pouvoir ordonner, le plus simple ici est de supprimer le premier character "J"

dsa_long_visit_numeric <-data.frame(dsa_long)
dsa_long_visit_numeric$VISIT<-as.numeric(substring(dsa_long_visit_numeric$VISIT,first = 2))


#on ordonne le tout:
dsa_long_sorted_by_numero_and_visit<-dsa_long_visit_numeric[order(dsa_long_visit_numeric$NUMERO,dsa_long_visit_numeric$VISIT),]
# on vérifie:
unique(dsa_long_sorted_by_numero_and_visit$NUMERO)
unique(dsa_long_sorted_by_numero_and_visit$VISIT)

# on vérifie aussi qu'on n'a plus de numéros dispérsés parmi les lignes:
for( i in unique(dsa_long_sorted_by_numero_and_visit[,1])){
  print(which(dsa_long_sorted_by_numero_and_visit[,1]==i))
}

# et qu'on a bien toutes nos observ:
dim(dsa_long_sorted_by_numero_and_visit)

dsa_long_sorted_by_numero_and_visit_rebootrowname<- dsa_long_sorted_by_numero_and_visit
rownames(dsa_long_sorted_by_numero_and_visit_rebootrowname)<-NULL
#d'abord on renome le df pour plus de lisibilité:
dsa_long_final<-dsa_long_sorted_by_numero_and_visit_rebootrowname

#Nous allons avoir recours à une liste, nous décidons de changer le numéro "0" pour "1"
dsa_long_final2<-dsa_long_final
dsa_long_final2[1,1]<-1


# Nous indexons les scores pour chaque numéros et chaques visites 

index_locf<-list()
for (i in unique(dsa_long_final2$NUMERO)){
  index_locf[[i]]<-unique(dsa_long_final2[dsa_long_final2$NUMERO==i,3])

}

#pour pouvoir ensuite calculer les sommes totales 




tot_locf <- vector(mode="list", length=length(index_locf))
for (i in 1:length(index_locf)){
for (j in 1:length(index_locf[[i]])){
 tot_locf[[i]][j] <-sum(dsa_long_final2[dsa_long_final2$NUMERO==i&dsa_long_final2$VISIT==index_locf[[i]][j],4:20])
}
}


#et prendre le score correspondant à la dernière visite. (On utilisera na.locf pour le NA à la dernière visite du numéro 128)

"%notin%"<-Negate("%in%")
library(zoo)

tot_locf_done_j56<-c()
for (i in 1:length(tot_locf)){
  if(i %notin% which(sapply(index_locf,is.null))){
  tot_locf_done_j56[i]<-tail(na.locf(unlist(tot_locf[i])),1)
  }
}
tot_locf_done_hdrs_j56<-na.omit(tot_locf_done_j56)
dsa_wide_sorted_by_numero<-dsa_wide[order(dsa_wide$NUMERO),]

# on remplace le numero 0 par 1:
dsa_wide_sorted_by_numero_0removed<-dsa_wide_sorted_by_numero
dsa_wide_sorted_by_numero_0removed[1,1]<-1

#on réinitialise là aussi les numéros de lignes:
dsa_wide_sorted_by_numero_0removed_rebootrowname<-dsa_wide_sorted_by_numero_0removed
rownames(dsa_wide_sorted_by_numero_0removed_rebootrowname)<-NULL
score_tot_hdrs_j0_sorted<-c()
for (i in unique(dsa_long_final2$NUMERO)){
score_tot_hdrs_j0_sorted[i]<-sum(dsa_long_final2[dsa_long_final2$NUMERO==i&dsa_long_final2$VISIT==0,4:20])
}
score_tot_hdrs_j0_sorted_naomit<-na.omit(score_tot_hdrs_j0_sorted)
dsa_wide_merged<-cbind(dsa_wide_sorted_by_numero_0removed_rebootrowname,score_tot_hdrs_j0_sorted_naomit,tot_locf_done_hdrs_j56)
check<-c()
for (i in 1:length(dsa_wide_sorted_by_numero_0removed_rebootrowname$NUMERO)){
  check[i]<-unique(dsa_wide_sorted_by_numero_0removed_rebootrowname$NUMERO)[i]==unique(dsa_long_final2$NUMERO)[i]
}

kable(table(check), caption="Match")
t.test(dsa_wide_merged[dsa_wide_merged[,2]==0,859],dsa_wide_merged[dsa_wide_merged[,2]==1,859],var.equal = T)
t.test(dsa_wide_merged[dsa_wide_merged[,2]==0,860],dsa_wide_merged[dsa_wide_merged[,2]==1,860],var.equal = T)
tot_par_visit<-c()
for (i in 1:nrow(dsa_long_final2)){
  tot_par_visit[i]<-sum(dsa_long_final2[i,4:20],na.rm = F)
}

length(tot_par_visit) # on a bien 1053 valeur qui correspond au nombre d'observations, dans l'ordre des lignes, en ayant gardé le NA

# on ajoute ce vecteur en tant que dernière colonne du Df:
dsa_long_final_tot_par_visit<-cbind(dsa_long_final2,tot_par_visit)

dsa_long_final_tot_par_visit_factor<-data.frame(dsa_long_final_tot_par_visit)

dsa_long_final_tot_par_visit_factor$VISIT<-as.numeric(as.factor(dsa_long_final_tot_par_visit_factor$VISIT))

#on se rassure: il n'ya pas eu de changement d'ordre dans cette manip, exemple: as.numeric(as.factor(c(1,98,32))) = [1] 1 3 2

library(lmerTest)

summary(lmer(formula = tot_par_visit~GROUPE+VISIT+(VISIT|NUMERO), data = dsa_long_final_tot_par_visit_factor))
library(ggplot2)
library(viridis)

ggplot(data=dsa_long_final_tot_par_visit_factor[dsa_long_final_tot_par_visit_factor$NUMERO %in%
                                                  unique(dsa_long_final_tot_par_visit_factor$NUMERO)[sample(x = unique(dsa_long_final_tot_par_visit_factor$NUMERO),size = 30,replace = F)],],
       mapping = aes(x = VISIT,y =tot_par_visit, group=as.factor(NUMERO),color=as.factor(NUMERO)) )+
  geom_line()+
  scale_colour_viridis(discrete = T)+
  theme(legend.position = "none")

ggplot(data =dsa_long_final_tot_par_visit_factor[dsa_long_final_tot_par_visit_factor$NUMERO%in%
                                                   unique(dsa_long_final_tot_par_visit_factor$NUMERO),],
       mapping = aes(x = VISIT,y = tot_par_visit,group=as.factor(NUMERO), color=as.factor(GROUPE)) )+
  geom_line()+
  scale_colour_viridis(discrete = T)+
  theme(legend.position = "none")
#Il nous faut créer un nouveau dataframne, pour ajouter un geom line spécifique

mean_group_visit<-aggregate(tot_par_visit~GROUPE+VISIT,dsa_long_final_tot_par_visit_factor,mean)


ggplot(data =dsa_long_final_tot_par_visit_factor[dsa_long_final_tot_par_visit_factor$NUMERO%in%
                                                   unique(dsa_long_final_tot_par_visit_factor$NUMERO),],
       mapping = aes(x = VISIT,y = tot_par_visit, color=as.factor(GROUPE)) )+
  geom_line(aes(group=as.factor(NUMERO)),alpha=0.4)+
  geom_line(data = mean_group_visit,size=3)+
  scale_colour_viridis(discrete = T)+
  theme(legend.position = "none")
cbinaire<-vector(mode="list", length=191)

for ( i in dsa_long_final_tot_par_visit_factor$NUMERO){
  for(j in dsa_long_final_tot_par_visit_factor[
    dsa_long_final_tot_par_visit_factor$NUMERO==i,3]){
    
      
   if(is.na(dsa_long_final_tot_par_visit_factor[
        dsa_long_final_tot_par_visit_factor$NUMERO==i&
          dsa_long_final_tot_par_visit_factor$VISIT==j,111])){
        cbinaire[[i]][j]<-NA
   }
    else if 
      (dsa_long_final_tot_par_visit_factor[
        dsa_long_final_tot_par_visit_factor$NUMERO==i&
        dsa_long_final_tot_par_visit_factor$VISIT==j,111] >
       dsa_long_final_tot_par_visit_factor[
         dsa_long_final_tot_par_visit_factor$NUMERO==i&
         dsa_long_final_tot_par_visit_factor$VISIT==1,111]/2){
        cbinaire[[i]][j]<-0
    }
     else{
       cbinaire[[i]][j]<-1
     }
  }
}


# Pour une raison indéterminée, 2 NA se sont glissés en plus dans les données, ce qui ramène les observations à 1055 au lieu de 1053. On les supprime:
which(sapply(X = cbinaire,FUN = function(X) NA %in% X))


cbinaire_f<-cbinaire
cbinaire_f[[92]]<-cbinaire_f[[92]][-2]
cbinaire_f[[170]]<-cbinaire_f[[170]][-3]


# Les résultats sont stockés en dernière colonne du data frame:

dsa_long_final_tot_par_visit_factor_2<-cbind(dsa_long_final_tot_par_visit_factor,unlist(cbinaire_f))
library(survival)

#Kaplan Meier
plot(survfit(Surv(VISIT,unlist(cbinaire_f))~GROUPE,data=dsa_long_final_tot_par_visit_factor_2),col=c("red","blue"))


#test du log-rank
survdiff(Surv(VISIT,unlist(cbinaire_f))~GROUPE,data=dsa_long_final_tot_par_visit_factor_2)