-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlm_across_windows.R
59 lines (52 loc) · 2.35 KB
/
lm_across_windows.R
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
## Linear model analyses across FC time windows
## FOR USE IN THE COGNITIVE AND BRAIN HEALTH LABORATORY
##################################################################################################################
##################################################################################################################
lm_across_window=function(window_duration,nthread, timeseries,model)
{
##check and remove subjects with incomplete data
idx.na=which(!complete.cases(model))
if(length(idx.na)>0)
{
model=model[-idx.na,]
timeseries=timeseries[-idx.na]
cat(paste0(length(idx.na)," subjects with missing data detected"))
}
##check number of frames
frames=rep(NA,length(timeseries))
for (sub in 1:length(timeseries)) {frames[sub]=NROW(timeseries[[sub]])}
nframes=min(frames)
cat(paste0("nframes=",nframes))
## correlation function to be used within lapply()
window.corr=function(ts,start,end)
{
cor.mat=cor(ts[start:end,])
cor.mat[upper.tri(cor.mat,diag = F)]
}
##activate parallel clusters
cl=parallel::makeCluster(nthread)
doParallel::registerDoParallel(nthread)
`%dopar%` = foreach::`%dopar%`
##looping across different window lengths
av.coef.vector=list()
for (window_length in 1:length(window_duration))
{
nwindows=nframes-window_duration[window_length]+1
av.coef.vector[[window_length]]=foreach::foreach(window=1:nwindows, .combine="c") %dopar%
{
#perform correlation analyses within each list element
all.cor.vector=lapply(timeseries, window.corr,start=1+window-1,end=window_duration[window_length]+window-1)
#flatten all list elements into a 2D matrix
all.cor.vector=scale(do.call(rbind,all.cor.vector))
mod=.lm.fit(y=all.cor.vector,x=data.matrix(cbind(1,model)))
coef=mod$coefficients[NCOL(model)+1,]
return(mean(abs(coef)))
}
}
return(av.coef.vector)
closeAllConnections()
}
##################################################################################################################
##################################################################################################################
#source("https://github.com/CogBrainHealthLab/MLtools/edit/main/lm_across_windows.R?raw=TRUE")
#lm_across_window(window_duration = c(80,160),nthread = 5,timeseries = CC.ts,model = cc.beh[,c("Sex","age_std")])