-
Notifications
You must be signed in to change notification settings - Fork 0
/
.Rhistory
194 lines (194 loc) · 7.46 KB
/
.Rhistory
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
??ggridges
library(ggplot2)
library(ggridges)
data <- data.frame(x = 1:5, y = rep(1, 5), height = c(0, 1, 3, 4, 2))
ggplot(data, aes(x, y, height = height)) + geom_ridgeline()
load("./ProjectsData.RData") # With this you dont need to run the next lines
load("~/gDrive/BCBL/PROYECTOS/2019 ThalamusNucleiStudy (Xabi)/ProjectsData.RData") # With this you dont need to run the next lines
GeneralPath ="~/gDrive/BCBL/PROYECTOS/2019 ThalamusNucleiStudy (Xabi)/results"
dir.create(GeneralPath)
GeneralPath ="/export/home/xintxausti/public/Gari/Xabi02/"
GeneralPath ="~/gDrive/BCBL/PROYECTOS/2019 ThalamusNucleiStudy (Xabi)/results"
dir.create(GeneralPath)
dimensions <- dim(ATLAS)
## Removing the Subjects with parkinson and dyslexia
eliminados <- 0
r <- c(0)
j=1;
for (i in 1:dimensions[1]) {
if (ATLAS[i,"Eliminar_1"] == 1 ){
r[j] <- c(i)
eliminados <- eliminados + 1
j=j+1
}
}
ATLAS <- ATLAS[-c(r),]
## Removing the Subject in which there is any movement in the MRI images (values 2 and 3)
remove <- 0
rr <- c(0)
j=1;
for (i in 1:dim(ATLAS)[1]) {
if (ATLAS[i, "DICOM_MOV"]!=1){
rr[j] <- c(i)
remove <- remove + 1
j=j+1
}
}
ATLAS <- ATLAS[-c(rr),]
## Removing the subject with bad segmentation
remove <- 0
rr <- c(0)
j=1;
for (i in 1:dim(ATLAS)[1]) {
if (ATLAS[i, "Bad_Segmentation"]==1){
rr[j] <- c(i)
remove <- remove + 1
j=j+1
}
}
ATLAS <- ATLAS[-c(rr),]
####
### NORMALIZATION
## Function to the normalization of the data, subnucleos with the correction extracting the effect of the ICV
myNorm <- function(x, y){
# x is the vector you want to normalize according to a linear regression
# over y, for example, x is hippocampal values and y is ICV (eTIV)
#??
#??Args.
# x: vector to normalize (hippocampal volume)
# y: vector to normalize it with (ICV)
# Returns:
# xn: vector x normalized
tempM <- na.omit(cbind(x, y))
xn <- x - lm(eval(x ~ y))$coefficients[2][[1]] * (y - mean(tempM[,2]))
return(xn)
}
## In this loop we normalize the dataset
ATLAS_Norm <- data.frame(ATLAS)
for (i in 1:52) {
ATLAS_Norm[,i] <- myNorm(ATLAS[,i], ATLAS$ICV)
}
ATLAS_Norm_Ord <- data.frame(matrix(0,dim(ATLAS_Norm)[1], dim(ATLAS_Norm)[2]))
colnames(ATLAS_Norm_Ord) <- t(nombresOrden)
for (i in 1:dim(ATLAS_Norm_Ord)[2]) {
for (j in 1:dim(ATLAS_Norm)[2]) {
if (colnames(ATLAS_Norm_Ord[i]) == colnames(ATLAS_Norm[j])) {
ATLAS_Norm_Ord[,i] = ATLAS_Norm[j]
}
}
}
##### Bootstrapping #####
# function to obtain R-Squared from the data
rsq <- function(formula, data, indices) {
d <- data[indices,] # allows boot to select sample
fit <- lm(formula, data=d, na.action = NULL)
#f = (100*(1-sum((d-fit$fitted.values)^2)/sum((d-mean(d))^2)));
return(summary(fit)$r.square)
#return(f)
}
bootstrapping <- function(data,n,PATH) {
##### LINEAR
# Initialise vectors to store statistics and CI limits in the for loop
# (actually it's not recommended to grow vectors inside loops but whatever)
nucNames <- colnames(data[1:n])
r2_linear <- vector()
lowlim_linear <- vector()
uplim_linear <- vector()
ResultadoBootsLinear <- list()
confintList_linear <- list()
predictions_linear <- list()
pvalor_lineal <- list()
malda <- list()
# Bootstrapping with 1000 replicates
for (i in 1:n) {
# bootstrapping, it is important to remove any NA, NaN or Inf values in the response variable
results <- boot(data = data[which(!is.na(data[,i])),], statistic=rsq,
R=1000, formula = as.formula(paste(nucNames[i], "AGE", sep="~")), sim="ordinary")
# get 95% confidence interval
confint <- boot.ci(results, type="norm")
# store everything for later
#ResultadoBootsLinear[[i]] <- results
ResultadoBootsLinear[[i]] <- results
confintList_linear[[i]] <- confint
r2_linear[i] <- confint$t0
lowlim_linear[i] <- confint$normal[2]
uplim_linear[i] <- confint$normal[3]
lineal_model <- lm(formula = as.formula(paste(nucNames[i], "AGE", sep="~")),data = data)
predictions_linear[[i]] <- lineal_model$fitted.values
pvalor_lineal[i] <- summary(lineal_model)$coefficients[2,4]
malda[i] <- confint$t0
}
###### Plot with 95% confidence intervals
#Red color for those not significant (CI95 includes 0)
colR2 <- lowlim_linear > 0
colR2[colR2 == "TRUE"] <- "blue"; colR2[colR2 == "FALSE"] <- "red"
#Plot with colors and custom axis
#jpeg(paste(PATH,"linear_R2.jpg", sep=""))
jpeg(height = 720, width = 1280, paste(PATH,"linear_R2.jpg", sep=""))
plotCI(1:n, r2_linear, li=lowlim_linear, ui=uplim_linear, xaxt="n", col=colR2,xlab="Nuclei", ylab="R2", pch=21, pt.bg=colR2)
abline(h=0, lty="dashed", col="gray")
text(1:n, par("usr")[3] - 0.01,
srt = 60, adj= 1, xpd = TRUE,
labels = nucNames[1:n], cex=0.65, col = colR2)
#legend(0, 0.58, legend=c("Significative", "Not significative"), col=c("blue", "red"), lty="solid", bty="n")
legend(45, 0.55, legend=c("Significative", "Not significative"), col=c("blue", "red"), lty="solid", bty="n")
dev.off()
# ## Se puede comenta
for (i in 1:n) {
title = paste(PATH,colnames(data[i]), sep="")
jpeg(height = 720, width = 1280, paste(title,"_Linear_Hist_R2.jpg", sep=""))
hist(ResultadoBootsLinear[[i]]$t, xlab="Bootstraped R2", main="Right LSg", col="lightblue", border="lightblue")
abline(v=ResultadoBootsLinear[[i]]$t0, lty="solid", lwd=2) #black vertical line, mean R2
abline(v=lowlim_linear[i], lty="dashed",lwd=2, col="red") #left red dashed line, lower limit of 95CI
abline(v=uplim_linear[i], lty="dashed",lwd=2, col="red") #right red dashed line, upper limit of 95CI
dev.off()
}
# for (i in 1:52) {
# title = paste(PATH,colnames(data[i]), sep="")
# jpeg(height = 720, width = 1280, paste(title,"_Linear_Density_R2.jpg", sep=""))
# plot(density(confint$t), lwd = 3, col = "steelblue", main = paste("Linear Model -",colnames(data[i])))
# abline(v = mean(confint$t), lwd=3, col = "gold")
# dev.off()
# }
###################
##### QUADRATIC ####
# Initialise vectors to store statistics and CI limits in the for loop
# (actually it's not recommended to grow vectors inside loops but whatever)
nucNames <- colnames(data[1:n])
r2_quadratic <- vector()
lowlim_quadratic <- vector()
uplim_quadratic <- vector()
ResultadoBootsQuadratic <- list()
confintList_quadratic <- list()
predictions_quadratic <- list()
pvalor_quadratic <- list()
# Bootstrapping with 1000 replicates
for (i in 1:n) {
# bootstrapping, it is important to remove any NA, NaN or Inf values in the response variable
results <- boot(data = data[which(!is.na(data[,i])),], statistic=rsq,
R=1000, formula = as.formula(paste(nucNames[i], "poly(AGE,2)", sep="~")), sim="ordinary")
# get 95% confidence interval
confint <- boot.ci(results, type="norm")
# store everything for later
ResultadoBootsQuadratic[[i]] <- results
confintList_quadratic[[i]] <- confint
r2_quadratic[i] <- confint$t0
lowlim_quadratic[i] <- confint$normal[2]
uplim_quadratic[i] <- confint$normal[3]
quadratic_model <- lm(formula = as.formula(paste(nucNames[i], "poly(AGE,2)", sep="~")),data = data)
predictions_quadratic[[i]] <- quadratic_model$fitted.values
pvalor_quadratic[i] <- summary(quadratic_model)$coefficients[3,4]
}
###### Plot with 95% confidence intervals
#Red color for those not significant (CI95 includes 0)
colR2 <- lowlim_quadratic > 0
colR2[colR2 == "TRUE"] <- "blue"; colR2[colR2 == "FALSE"] <- "red"
#Plot with colors and custom axis
jpeg(height = 720, width = 1280, paste(PATH,"quadratic_R2.jpg", sep=""))
plotCI(1:n, r2_quadratic, li=lowlim_quadratic, ui=uplim_quadratic, xaxt="n", col=colR2,xlab="Nuclei", ylab="R2", pch=21, pt.bg=colR2)
abline(h=0, lty="dashed", col="gray")
text(1:n, par("usr")[3] - 0.01,
srt = 60, adj= 1, xpd = TRUE,
labels = nucNames[1:n], cex=0.65, col = colR2)
legend(45, 0.55, legend=c("Significative", "Not significative"), col=c("blue", "red"), lty="solid", bty="n")
dev.off()