-
Notifications
You must be signed in to change notification settings - Fork 0
/
LongitInclude.cpp
1180 lines (1016 loc) · 43.1 KB
/
LongitInclude.cpp
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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include <iostream>
#include <string>
#include <RcppArmadillo.h>
#include <Rcpp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace arma;
using namespace Rcpp;
using namespace std;
////////////////////////////// Functions used in cross sectional and Longitudinal Models //////////
//calculate cross covariance for an individual visit
// [[Rcpp::export]]
cube ccv(const mat& xt, int lag){
//xt is matrix for a given visit moving across ROIs and down time
//lag is maximum number of lags to use in the calculations
int P = xt.n_cols; //number of ROIs
int T = xt.n_rows; //length of time series
cube covMat(P, P, 2*lag+1, fill::zeros);
vec xi(T, fill::zeros);
double xibar;
vec xj(T, fill::zeros);
double xjbar;
for(int i = 0; i < P; i++){
for(int j = 0; j < P; j++){
xi = xt.col(i);
xj = xt.col(j);
xibar=mean(xi);
xjbar=mean(xj);
for(int l = 0; l <= lag; l++){
vec xil = xi.head(T-l);
vec xjl = xj.tail(T-l);
double sum = 0;
for(int t=0; t<(T-l); t++){
sum += (xil(t)-xibar)*(xjl(t)-xjbar);
}
covMat(i,j,lag+l) = covMat(j,i,lag-l) = sum/T;//cross covariance for ROI i and j at lag l
}
}
}
return covMat; //return cube with ROIs in rows and columns and lags in slices (lag 0 is middle slice)
}
//function to calculate the lag 0 cross correlation matrix for a given visit
// [[Rcpp::export]]
mat ccr(const cube& C){
//C is a PxPxL cross covariance cube output by the ccv function
int P = C.n_cols; //Number of ROIs
int L = C.n_slices; //2 times the number of lags plus 1
mat R(P,P,fill::zeros);
mat lag0cov = C.slice((L-1)/2); //take the lag 0 slice from covariance cube
for(int i = 0; i < P; i++){
for(int j = 0; j < P; j++){
R(i,j)=lag0cov(i,j)/std::sqrt(lag0cov(i,i)*lag0cov(j,j)); //convert to correlation
}
}
return R; //return correlation matrix for lag 0
}
//calculate weights to be used in the estimation of delta
// [[Rcpp::export]]
arma::vec MB(double lag, double b){
//lag is maximum number of lags to use in calculation
//b is the bandwidth of the windowing function times sqrt(T) where T is the length of the time series
vec weights(2*lag+1, fill::zeros);
for(int i = -lag; i < lag+1; i++){
if(std::abs(i/b)<=1){weights[i+lag] = 1-std::abs(i/b);}
else{weights[i+lag] = 0;}
}
return weights; //return vector of weights of length 2 times lag + 1 (lag 0 in center)
}
//function to calculate delta which are used in calculation of Roy covariance cov(r_ij, r_lm)
// [[Rcpp::export]]
double delta(int i, int j, int l, int m, const vec& weights, const cube& C, int lag, int lagln){
//i, j, l, and m are the ROIs for the desired delta value
//weights is a vector of output weights from the MB function
//C is cross covariance cube from ccv function
//lag is maximum lag to be used in calculations
//lagln is 2*lag+1
double sum=0;
for(int t=0; t<lagln; t++){
sum += weights(t)*C(i-1,j-1,t)*C(l-1,m-1,t); //estimate theta
}
double result = sum/std::sqrt(C(i-1,i-1,lag)*C(j-1,j-1,lag)*C(l-1,l-1,lag)*C(m-1,m-1,lag)); //convert to correlation
return result;
}
//function calculate Roy covariance cov(r_ij, r_lm)
// [[Rcpp::export]]
double roycov(int i, int j, int l, int m, const cube& C, const mat& R, double T,
const vec& weights, int lag, int lagln){
//i,j,l,m are ROIs of the desired covariance
//C is cross covariance cube from ccv function
//R is lag 0 cross correlation from ccr function
//T is the length of the time series
//weights is a vector of weights from the MB function(or other windowing function)
//lag is maximum lag to be used in calculations
//lagln is 2*lag+1
//formula for Roy covariance
double result=(0.5*R(i-1,j-1)*R(l-1,m-1)*(delta(i,l,i,l,weights,C,lag,lagln)+
delta(i,m,i,m,weights,C,lag,lagln)+
delta(j,l,j,l,weights,C,lag,lagln)+
delta(j,m,j,m,weights,C,lag,lagln)) -
R(i-1,j-1)*(delta(i,l,i,m,weights,C,lag,lagln)+delta(j,l,j,m,weights,C,lag,lagln)) -
R(l-1,m-1)*(delta(j,l,i,l,weights,C,lag,lagln)+delta(j,m,i,m,weights,C,lag,lagln)) +
delta(i,l,j,m,weights,C,lag,lagln) + delta(j,l,i,m,weights,C,lag,lagln))/T;
return result;
}
//construct within subject covariance matrix using Roy covariance function, unstructured
// [[Rcpp::export]]
mat RoyMat(const cube& C, const mat& R, double T, const vec& weights, int lag, string SigmaType){
//C is cross covariance cube from ccv function
//R is lag 0 cross correlation matrix from ccr function
//T is the length of the time series
//weights is a vector of weights from the MB function(or other windowing function)
//lag is maximum lag to be used in calculations
//SigmaType is a string to signify if SigmaRoy variance should be "Zero", "Unstructured", or "Diagonal"
int P=C.n_rows; //number of ROIs
int Q=P*(P-1)/2; //number of ROI pairs
int lagln=2*lag+1;
mat result(Q,Q,fill::zeros);
if(SigmaType=="Unstructured"){
vec resultvec(Q*Q, fill::zeros);
//create vector of roy covariances
int count=0;
for(int i=0; i<(P-1); i++){
for(int j=i+1; j<P; j++){
for(int l=0; l<(P-1); l++){
for(int m=l+1; m<P; m++){
if((l > i) || (l == i && m >= j)){
resultvec(count) = roycov(i+1,j+1,l+1,m+1,C,R,T,weights,lag,lagln);
count += 1;
}
}
}
}
}
//populate within subject Roy variance matrix from the vector of Roy variances
for(int i=0; i<Q; i++){
for(int j=i; j<Q; j++){
result(i,j)=resultvec[0];
result(j,i)=resultvec[0];
resultvec=resultvec.tail(resultvec.n_elem-1);
}
}
}
if(SigmaType=="Diagonal"){
int count=0;
for(int i=0; i<(P-1); i++){
for(int j=i+1; j<P; j++){
result(count,count) = roycov(i+1,j+1,i+1,j+1,C,R,T,weights,lag,lagln);
count += 1;
}
}
}
if(SigmaType=="Zero"){
result.fill(0);
}
return result;
}
//Apply previous functions to get within subject Roy variance for a single visit
// [[Rcpp::export]]
List Roy(const mat& x, int lag, double bw, string SigmaType){
//x is matrix for a given subject moving across ROIs and down time
//lag is maximum number of lags to use in the calculations
//SigmaType is a string to signify if SigmaRoy variance should be "Zero", "Unstructured", or "Diagonal"
//bw is the bandwidth to be used by the MB function
int T=x.n_rows; //length of time series
lag=min(lag, T);
cube C=ccv(x,lag); //cross covariance cube
mat R=ccr(C); //lag 0 correlation matrix
double b = bw*sqrt(T); //bandwidth used in windowing function
vec weights=pow(MB(lag,b),2);
List out(2);
out("Sigma")=RoyMat(C,R,T,weights,lag,SigmaType);
out("R")=R;
return out;
}
//function to calculate test statistic for vector beta
// [[Rcpp::export]]
double getStat(const vec& beta1, const vec& beta2, const mat& beta1var, const mat& beta2var){
//beta1 and beta1var are the estimates and variances of beta for group 1
//beta2 and beta2var are the estimates and variances of beta for group 2
mat Contrast=eye(beta1.n_elem, beta1.n_elem);
mat VarComb = beta1var + beta2var;
vec betadiff = beta1 - beta2;
mat T(1,1);
//calculate test statistic
T = trans(Contrast*betadiff) * inv(Contrast*VarComb*Contrast) * Contrast*betadiff;
return T(0,0);
}
//function to calculate test statistic for scaler beta
// [[Rcpp::export]]
double getStatUni(double beta1, double beta2, double beta1var, double beta2var){
//beta1 and beta1var are the estimates and variances of beta for group 1
//beta2 and beta2var are the estimates and variances of beta for group 2
return pow(beta1 - beta2, 2) / (beta1var + beta2var);
}
//function to estimate Psi0
// [[Rcpp::export]]
mat get_Psi0(const mat& error, const cube& RoyVars, int Visits, int Q, string Psi0type){
//error is a residual matrix where each column represents a visit
//RoyVars is a cube with each slice being the Roy within visit covariance matrix for one visit
//Visits is the total number of visits (or number of subjects in cross sectional model)
//Q is the number of ROI pairs
//Psi0type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
mat Empir = error*error.t()/Visits; //empirical error covariance matrix given last beta
mat AvgRoy = mean(RoyVars,2); //average Roy variance matrix across all N subjects
mat Psi0(Q, Q);
if(Psi0type == "Zero"){//Solve for Psi0 if it is assumed to be 0
Psi0.fill(0);
}
if(Psi0type == "Unstructured"){//Solve for Psi0 if unstructured structure is assumed
for(int i = 0; i < Q; i++){
for(int j = 0; j < Q; j++){
Psi0(i,j) = Empir(i,j) - AvgRoy(i,j);
}
}
for(int i = 0; i < Q; i++){
if(Psi0(i,i) <= 0){Psi0(i,i) = 0;}
}
}
if(Psi0type == "Diagonal"){//Solve for Psi0 if diagonal structure is assumed
for(int i = 0; i < Q; i++){
for(int j = 0; j < Q; j++){
if(i == j){Psi0(i,j) = Empir(i,j) - AvgRoy(i,j);}else{Psi0(i,j) = 0;}
}
}
for(int i = 0; i < Q; i++){
if(Psi0(i,i) <= 0){Psi0(i,i) = 0;}
}
}
if(Psi0type == "Scaled"){//Solve for Psi0 if scaled identity structure is assumed
vec diffsigma2(Q);
double sigma2;
for(int i = 0; i < Q; i++){
diffsigma2(i) = Empir(i,i) - AvgRoy(i,i);
}
if(mean(diffsigma2) < 0){sigma2 = 0;}else{sigma2 = mean(diffsigma2);}
for(int i = 0; i < Q; i++){
for(int j = 0; j < Q; j++){
if(i == j){Psi0(i,j) = sigma2;}else{Psi0(i,j) = 0;}
}
}
}
if(Psi0type == "CS"){//Solve for Psi0 if compound symmetry structure is assumed
double kappa=0;
if(Q != 1){
vec diffkappa(Q*(Q-1)/2);
int count = 0;
for(int i = 0; i < Q; i++){
for(int j = i+1; j < Q; j++){
diffkappa(count) = Empir(i,j) - AvgRoy(i,j);
count++;
}
}
kappa = mean(diffkappa);
}
vec diffsigma2(Q);
double sigma2;
for(int i = 0; i < Q; i++){
diffsigma2(i) = Empir(i,i) - AvgRoy(i,i);
}
if(mean(diffsigma2) < 0){sigma2 = 0;}else{sigma2 = mean(diffsigma2);}
for(int i = 0; i < Q; i++){
for(int j = 0; j < Q; j++){
if(i == j){Psi0(i,j) = sigma2;}else{Psi0(i,j) = kappa;}
}
}
}
return Psi0;
}
///////////////////////// Functions used in cross sectional model //////////////////
//function to estimate beta assuming diagonal design matrix in cross sectional model
// [[Rcpp::export]]
List getBeta(int N, int Q, const cube& invTotVar, const mat& r){
//N is number of subjects in the group
//Q is number of ROI pairs
//invTotVar is a cube where each slice is the inverse of Roy+Psi for one subject
//r is an NxQ matrix of the observed correlations
mat SumInvVar=sum(invTotVar,2); //sum up inverse covariance matrices
mat betaVar=inv(SumInvVar); //invert to estimate the variance matrix of beta
//estimate beta
cube temp(Q,1,N,fill::zeros);
for(int n = 0; n < N; n++){
temp.slice(n)=invTotVar.slice(n)*r.row(n).t();
}
vec tempsum=sum(temp,2);
vec beta=betaVar*tempsum;
return List::create(
_("betaVar")=betaVar,
_("beta")=beta
);
}
//function to iteratively update beta and Psi using method of moments for cross sectional model
// [[Rcpp::export]]
List getPsiBeta(int N, int Q, const mat& r, const cube& RoyVars, double step2_tol,
string Psi0type, int MaxIter){
//N is number of subjects in the group
//Q is number of ROI pairs
//r is matrix of observed correlations with a row for each subject
//RoyVars is a cube where each slice is the estimated Roy within subject variance for one subject
//step2_tol is used for the stopping criteria
//Psi0type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//MaxIter is maximum iterations to perform
bool done = FALSE;
int iter = 0; // Current iteration
int code = 0; // 0 - converged, 1 - otherwise
//initiate estimates
List betaout(2);
vec beta(Q,fill::zeros);
vec betaold(Q,fill::zeros);
mat Psi0(Q,Q,fill::zeros);
mat Psi0old(Q, Q, fill::zeros);
cube TotVar(Q,Q,N,fill::zeros);
cube invTotVar(Q,Q,N,fill::zeros);
betaold=mean(r,0).t();
mat error(Q, N, fill::zeros);
for(int n = 0; n < N; n++){
for(int q = 0; q < Q; q++){
error(q,n) = r(n,q) - betaold(q);
}
}
//iterate until convergence or max number of iterations hit
while(!done){
for(int n = 0; n < N; n++){
for(int i = 0; i < Q; i++){
error(i,n) = r(n,i) - beta(i);
}
}
Psi0=get_Psi0(error, RoyVars, N, Q, Psi0type=Psi0type); //update Psi0
for(int n = 0; n < N; n++){
TotVar.slice(n)=RoyVars.slice(n)+Psi0;
invTotVar.slice(n)=inv(TotVar.slice(n));
}
betaout=getBeta(N, Q, invTotVar, r);
beta=as<vec>(betaout("beta")); //update beta
for(int n = 0; n < N; n++){
for(int q = 0; q < Q; q++){
error(q,n) = r(n,q) - beta(q);
}
}
//stop if change in estimates is small
if((norm(beta-betaold,2) < step2_tol) &&
(abs(Psi0 - Psi0old).max() < step2_tol)){
done = TRUE;
code = 0;
}
iter++;
//stop if max number of iterations is reached
if(iter >= MaxIter){
done = TRUE;
code = 1;
}
betaold = beta;
Psi0old = Psi0;
}
return List::create(
_["beta"] = beta,
_["betaVar"] = as<mat>(betaout("betaVar")),
_["Psi0"] = Psi0,
_["TotVar"] = TotVar,
_["code"] = code
);
}
//function to permute data for permutation test for cross sectional model
// [[Rcpp::export]]
List permute(int N, int Q, const mat& r, const cube& RoyVars, vec ids){
//N is number of subjects in both groups combined
//Q is number of ROI pairs
//r is observed correlations with a row for each subject
//RoyVars is a cube of within subject Roy variances with a slice for each subject
//ids is a vector containing (0, 1,... , n-1)
vec ids_new=shuffle(ids); //permute indices
mat r_new(N,Q);
cube RoyVars_new(Q,Q,N);
int id;
for(int n = 0; n < N; n++){
id=ids_new(n);
r_new.row(n)=r.row(id); //reorder r
RoyVars_new.slice(n)=RoyVars.slice(id); //reorder RoyVars
}
return List::create(
_["RoyVars_new"] = RoyVars_new,
_["r_new"] = r_new
);
}
//function to permute data and calculate test statistic for permuted data for cross sectional model
// [[Rcpp::export]]
List permT(int N1, int N2, int Q, const mat& r, const cube& RoyVars, double step2_tol,
vec ids, string Psi0type, bool verbose, int MaxIter){
//N1 is number of subjects in group 1
//N2 is number of subjects in group 2
//Q is number of ROI pairs
//r is observed correlations with a row for each subject
//RoyVars is a cube of within subject Roy variances with a slice for each subject
//step2_tol is used for stoping criteria for solving for beta and Psi
//ids is a vector containing (0, 1,... , n-1)
//Psi0type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//verbose is a boolean to denote if permutation progress updates should be printed
//MaxIter is the maximum number of IWLS iterations to run
int N=N1+N2;
List permdat=permute(N, Q, r, RoyVars, ids); //permute the data
mat r_new=permdat("r_new");
cube RoyVars_new=permdat("RoyVars_new");
//recalculate beta1 and beta1var with new group assignment
List PsiBeta1_new=getPsiBeta(N1, Q, r_new.head_rows(N1), RoyVars_new.head_slices(N1),
step2_tol, Psi0type, MaxIter);
vec beta1_new=PsiBeta1_new("beta");
mat beta1Var_new=PsiBeta1_new("betaVar");
//recalculate beta2 and beta2var with new group assignment
List PsiBeta2_new=getPsiBeta(N2, Q, r_new.tail_rows(N2), RoyVars_new.tail_slices(N2),
step2_tol, Psi0type, MaxIter);
vec beta2_new=PsiBeta2_new("beta");
mat beta2Var_new=PsiBeta2_new("betaVar");
//recalculate test statistic with new group assignments
double T_global=getStat(beta1_new, beta2_new, beta1Var_new, beta2Var_new);
vec T_local(Q, fill::zeros);
for(int q = 0; q < Q; q++){
T_local(q) = getStatUni(beta1_new(q), beta2_new(q),
beta1Var_new(q, q), beta2Var_new(q, q));
}
return List::create(
_["T_global"] = T_global,
_["T_local"] = T_local
);
}
//master function for cross sectional model
// [[Rcpp::export]]
List FCanalysis(const List& datalist, int N1, int N2, int Nperms, int lag=50,
double bw=5, double step2_tol=1e-5, string Psi0type="CS", int MaxIter=20,
bool verbose=false, string SigmaType="Unstructured") {
//datalist is a list where each element is a TxP matrix of data for one subject
//N1 is the number of subjects in group 1
//N2 is the number of subjects in group 2
//Nperms is the number of permutations to run for the permutation test
//lag is the maximum number of lags to be used in calculations
//bw is the bandwidth to be used for the MB windowing function
//step2_tol is the stoping criteria for iteratively solve for beta and Psi
//Psi0type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//MaxIter is maximum iterations to perform solving for Psi and beta
//verbose is boolean to signify if permutation test status updates should be printed to the console
//SigmaType is a string to signify if SigmaRoy variance should be "Zero", "Unstructured", or "Diagonal"
///Define dimension parameters and output list
List outList;
int N = N1+N2;
int P=as<mat>(datalist(0)).n_cols;
int Q=P*(P-1)/2;
///Estimate Roy Variances
cube RoyVars(Q,Q,N);
cube R(P,P,N);
for(int n = 0; n < N; n++) {
List Royn=Roy(datalist(n), lag, bw, SigmaType);
RoyVars.slice(n) = as<mat>(Royn("Sigma"));
R.slice(n) = as<mat>(Royn("R"));
}
cube RoyVarsG1(Q,Q,N1);
RoyVarsG1=RoyVars.head_slices(N1);
cube RoyVarsG2(Q,Q,N2);
RoyVarsG2=RoyVars.tail_slices(N2);
outList("SigmaG1")=RoyVarsG1;
outList("SigmaG2")=RoyVarsG2;
////build r estimate matrix
mat r(N,Q);
for(int n = 0; n < N; n++){
int count=0;
for(int j = 0; j < P; j++){
for(int k = j+1; k < P; k++){
r(n,count)=R(j,k,n);
count++;
}
}
}
mat rG1=r.head_rows(N1);
mat rG2=r.tail_rows(N2);
outList("rG1")=rG1;
outList("rG2")=rG2;
///Estimate Psi and beta for both groups
List PsiBetaG1=getPsiBeta(N1, Q, rG1, RoyVarsG1, step2_tol, Psi0type, MaxIter);
vec betaG1=PsiBetaG1("beta");
mat betaG1Var=PsiBetaG1("betaVar");
mat Psi0G1=PsiBetaG1("Psi0");
cube TotVarG1=PsiBetaG1("TotVar");
double codeG1=PsiBetaG1("code");
if((codeG1==1) & (MaxIter>1)){cout << "WARNING: Group 1 estimates did not converge in " << MaxIter << " iterations" << endl;}
List PsiBetaG2=getPsiBeta(N2, Q, rG2, RoyVarsG2, step2_tol, Psi0type, MaxIter);
vec betaG2=PsiBetaG2("beta");
mat betaG2Var=PsiBetaG2("betaVar");
mat Psi0G2=PsiBetaG2("Psi0");
cube TotVarG2=PsiBetaG2("TotVar");
double codeG2=PsiBetaG2("code");
if((codeG2==1) & (MaxIter>1)){cout << "WARNING: Group 2 estimates did not converge in " << MaxIter << " iterations" << endl;}
outList("Psi0G1")=Psi0G1;
outList("Psi0G2")=Psi0G2;
outList("TotalVarianceG1")=TotVarG1;
outList("TotalVarianceG2")=TotVarG2;
outList("betaG1var")=betaG1Var;
outList("betaG1")=betaG1;
outList("betaG2var")=betaG2Var;
outList("betaG2")=betaG2;
//Calculate Test Statistics
double T_global;
T_global=getStat(betaG1, betaG2, betaG1Var, betaG2Var);
outList("T_global")=T_global;
vec T_local(Q);
for(int q = 0; q < Q; q++){
T_local(q) = getStatUni(betaG1(q), betaG2(q), betaG1Var(q, q), betaG2Var(q, q));
}
outList("T_local")=T_local;
//Run permutation test
vec T_dist_global(Nperms, fill::zeros);
mat T_dist_local(Nperms, Q);
vec ids(N,fill::zeros);
for(int n = 0;n < N; n++){
ids(n)=n;
}
List out2(2);
for(int nperm = 0; nperm < Nperms; nperm++){
out2 = permT(N1, N2, Q, r, RoyVars, step2_tol, ids, Psi0type, verbose, MaxIter);
T_dist_global(nperm) = as<double>(out2("T_global"));
T_dist_local.row(nperm) = as<vec>(out2("T_local")).t();
if(verbose==true){cout << "Permutation " << nperm+1 << " of " << Nperms << endl;}
}
outList("T_dist_global")=T_dist_global;
outList("T_dist_local")=T_dist_local;
return outList;
}
///////////////////////////////// Functions used in Longitudinal Models //////////////////////////
//Function to estimate Psi1 for longitudinal model
// [[Rcpp::export]]
mat get_Psi1(const vec& visvec, int Visits, int maxVis, int N, int Q, const mat& error,
string Psi1type){
//visvec is a vector with the number of visits for each subject
//Visits is the number of total visits, i.e. sum(visvec)
//maxVis is the maximum number of visits by any one subject
//N is the number of subjects
//Q is the number of ROI pairs
//error is a residual matrix where each column represents a visit
//Psi1type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//create a cube of Psi1 estimates from each possible combination or visits
cube Psi1cube(Q, Q, maxVis*(maxVis-1)/2, fill::zeros);
int countouter = 0;
int countinner;
mat v1mat;
mat v2mat;
vec v2ind(N, fill::zeros);
int Nv2Tot = 0;
int Nv2;
mat AvgPsi1(Q, Q, fill::zeros);
if(Psi1type != "Zero"){
for(int v1 = 0; v1 < maxVis-1; v1++){
for(int v2 = v1+1; v2 < maxVis; v2++){
for(int n = 0; n < N; n++){
if(visvec(n) > v2){v2ind(n) = 1;}else{v2ind(n) = 0;}
}
Nv2 = sum(v2ind); //number of subjects with a visits at v1 and v2
Nv2Tot += Nv2;
v1mat.resize(Q, Nv2);
v2mat.resize(Q, Nv2);
countinner = 0;
for(int n = 0; n < N; n++){
if(v2ind(n) == 1){
v1mat.col(countinner) = error.col(sum(visvec.head(n))+v1);
v2mat.col(countinner) = error.col(sum(visvec.head(n))+v2);
countinner++;
}
}
Psi1cube.slice(countouter) = v1mat*v2mat.t();
countouter++;
}
}
AvgPsi1=sum(Psi1cube,2)/Nv2Tot;
}
mat Psi1(Q,Q,fill::zeros); //Initialize Psi1
if(Psi1type == "Zero"){
Psi1.fill(0);
}
if(Psi1type == "CS"){ //Solve for Psi1 if compound symmetry structure for Psi1 is assumed
double omega; //Initialize omega, the off diagonal element of Psi1
double nu; //Initialize nu, the diagonal element of Psi1
if(Q == 1){omega = 0;} //if P=2 then Q=1 and Psi1 is a scaler so there are no off diagonal elements
else{
vec omegavec(Q*(Q-1)/2);
int count = 0;
for(int i = 0; i < Q; i++){
for(int j = i+1; j < Q; j++){
omegavec(count) = AvgPsi1(i,j);
count++;
}
}
omega = mean(omegavec); //omega is mean of all off diagonal values of all slices of Psi1cube
}
vec nuvec(Q);
for(int i = 0; i < Q; i++){
nuvec(i) = AvgPsi1(i,i);
}
nu = mean(nuvec); //nu is mean of all diagonal values of all slices of Psi1cube
for(int i = 0; i < Q; i++){
for(int j = 0; j < Q; j++){
if(i == j){Psi1(i,j) = nu;}else{Psi1(i,j) = omega;}
}
}
}
if(Psi1type == "Scaled"){ //Solve for Psi1 if scaled identity structure for Psi1 is assumed
double nu; //Initialize nu, the diagonal element of Psi1
vec nuvec(Q);
for(int i = 0; i < Q; i++){
nuvec(i) = AvgPsi1(i,i);
}
nu = mean(nuvec); //nu is mean of all diagonal values of all slices of Psi1cube
for(int i = 0; i < Q; i++){
Psi1(i,i) = nu;
}
}
if(Psi1type == "Unstructured"){//Solve for Psi1 if unstructured structure is assumed
Psi1=AvgPsi1;
}
if(Psi1type == "Diagonal"){//Solve for Psi1 if diagonal structure is assumed
Psi1=AvgPsi1;
for(int i = 0; i < Q; i++){
for(int j = 0; j < Q; j++){
if(i != j){Psi1(i,j) = 0;}
}
}
}
return Psi1;
}
//function to iteratively update beta and Psi using method of moments for longitudinal model
// [[Rcpp::export]]
List getPsiBetaLongit(int N, int Q, const mat& r, const cube& RoyVars, double step2_tol,
int MaxIter, const vec& visvec, const mat& X,
string Psi0type, string Psi1type, int Visits, int maxVis){
//N is number of subjects in the group
//Q is number of ROI pairs
//r is matrix of observed correlations with a row for each visit
//RoyVars is a cube where each slice is the estimated Roy within visit variance for one visit
//step2_tol is used for the stopping criteria
//MaxIter is maximum iterations to perform
//visvec is a vector with the number of visits for each subject
//X is the design matrix
//Psi0type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//Psi1type denotes the structure to be used for Psi0, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//Visits is the total number of visits, i.e. sum visvec
//maxVis is the maximum number of visits for any one subject
bool done = FALSE;
int iter = 0; // Current iteration
int code = 0; // 0 - converged, 1 - otherwise
//Initiate beta estimates at the OLS estimates
vec Y = vectorise(r.t());
vec betaold = inv(X.t() * X) * X.t() * Y;
mat fits = X * betaold;
mat error(Q, Visits, fill::zeros);
for(int v = 0; v < Visits; v++){
for(int q = 0; q < Q; q++){
error(q,v) = r(v,q) - fits(v*Q+q);
}
}
//Initiate other parameters
vec beta(betaold.n_elem, fill::zeros);
mat betaVar(betaold.n_elem, betaold.n_elem, fill::zeros);
mat Psi0old(Q, Q, fill::zeros);
mat Psi1old(Q, Q, fill::zeros);
mat Psi0(Q, Q);
mat Psi1(Q, Q);
mat TotVar(Q*Visits, Q*Visits, fill::zeros);
mat InvTotVar(Q*Visits, Q*Visits, fill::zeros);
//iterate until convergence or max number of iterations hit
while(!done){
Psi0=get_Psi0(error, RoyVars, Visits, Q, Psi0type); //update Psi0
Psi1=get_Psi1(visvec, Visits, maxVis, N, Q, error, Psi1type); //update Psi1
int counter=0;
for(int n = 0; n < N; n++){
for(int v = 0; v < visvec(n); v++){
TotVar.submat((counter+v)*Q, (counter+v)*Q, size(Q,Q))=RoyVars.slice(counter+v)+Psi0;
for(int v2 = v+1; v2 < visvec(n); v2++){
TotVar.submat((counter+v)*Q, (counter+v2)*Q, size(Q,Q)) = Psi1;
TotVar.submat((counter+v2)*Q, (counter+v)*Q, size(Q,Q)) = Psi1;
}
}
InvTotVar.submat(counter*Q, counter*Q, size(Q*visvec(n), Q*visvec(n))) =
inv(TotVar.submat(counter*Q, counter*Q, size(Q*visvec(n), Q*visvec(n))));
counter += visvec(n);
}
betaVar = inv(X.t()*InvTotVar*X); //Estimate variance of beta
beta = betaVar*X.t()*InvTotVar*Y; //Update beta
fits = X * beta;
for(int v = 0; v < Visits; v++){
for(int q = 0; q < Q; q++){
error(q,v) = r(v,q) - fits(v*Q+q);
}
}
//stop if change in estimates is small
if((norm(beta-betaold,2) < step2_tol) &&
(abs(Psi0 - Psi0old).max() < step2_tol) &&
(abs(Psi1 - Psi1old).max() < step2_tol)){
done = TRUE;
code = 0;
}
iter++;
//stop if max number of iterations is reached
if(iter >= MaxIter){
done = TRUE;
code = 1;
}
betaold = beta;
Psi0old = Psi0;
Psi1old = Psi1;
}
return List::create(
_["beta"] = beta,
_["betaVar"] = betaVar,
_["Psi0"] = Psi0,
_["Psi1"] = Psi1,
_["TotVar"] = TotVar,
_["code"] = code,
_["Residuals"] = error.t()
);
}
//function to resample the data for bootstrap or permutation test
// [[Rcpp::export]]
List Resample(int N, int Q, const mat& r, const cube& RoyVars, const mat& Xall,
const vec& visvec, bool replace, const vec& ids){
//N is the number of subjects
//Q is the number of ROI pairs
//r is a matrix of observed correlation coefficients with a row for each visit
//RoyVars is a cube of estimated within visit Roy variances with a slice for each visit
//Xall is the design matrices from both groups stacked
//visvec is a vector with the number of visits for each subject
//replace is a boolean, if false then data is permuted, if true then data is bootstraped
//ids is a vector with integer values 0 through N-1
vec ids_new(N);
if(replace == true){
vec randoms(N);
randoms.randu();
ids_new = floor(randoms*N);
}
if(replace == false){
ids_new = shuffle(ids); //permute indices
}
vec visvec_new(N);
for(int n = 0; n < N; n++){
visvec_new(n) = visvec(ids_new(n));
}
int SampVisits=sum(visvec_new);
double counter=0;
double lookup=0;
mat r_new(SampVisits,Q);
cube RoyVars_new(Q,Q,SampVisits);
mat Xall_new(SampVisits*Q, Xall.n_cols);
for(int n = 0; n < N; n++){
lookup = sum(visvec.head(ids_new(n)+1)) - visvec(ids_new(n));
for(int v = 0; v < visvec_new(n); v++){
r_new.row(counter+v) = r.row(lookup+v); //reorder r
RoyVars_new.slice(counter+v) = RoyVars.slice(lookup+v); //reorder RoyVars
Xall_new.submat((counter+v)*Q, 0, size(Q,Xall.n_cols)) = Xall.submat((lookup+v)*Q, 0, size(Q,Xall.n_cols));
}
counter+=visvec_new(n);
}
return List::create(
_["RoyVars_new"] = RoyVars_new,
_["r_new"] = r_new,
_["Xall_new"] = Xall_new,
_["visvec_new"] = visvec_new
);
}
//function to calculate test statistic on permuted full model residuals for multivariate test
// [[Rcpp::export]]
List PermTestLongit(int N1, int N2, int Q, const mat& r_resid_new, const cube& RoyVars_new,
const mat& Xall_new, const vec& visvec_new, double step2_tol,
string Psi0type, string Psi1type, int MaxIter,
const vec& betaG1, const vec& betaG2, int start, int end,
int Group1Visits_new, int Group2Visits_new, int Group1maxVis_new,
int Group2maxVis_new){
//N1 and N2 are the number of subjects in group 1 and group 2 respectively
//Q is number of ROI pairs
//r_resid_new are permuted full model residuals with a row for each visit
//RoyVars_new is permuted a cube of within visit Roy variances with a slice for each visit
//Xall_new is the permuted design matrices from the two groups stacked, output from Resample
//visvec_new is a permuted vector with the number of visits for each subject
//step2_tol is used for stoping criteria for solving for beta and Psi
//Psi0type and Psi1type denote the structures to be used for Psi0 and Psi1, either "CS" for compound symmetry,
//"Scaled" for Scaled Identity, "Zero", "Unstructured", or "Diagonal"
//MaxIter is the maximum number of IWLS iterations to run
//betaG1 and betaG2 are vectors of estimated coefficients from full model for groups 1 and 2 respectively
//start and end are the first and last indices of beta involved in the test
//Group1Visits_new and Group2Visits_new are the number of visits in the permuted data for groups 1 and 2 respectively
//Group1maxVis_new and Group2maxVis_new are the maximum number of visits in the permuted data for groups 1 and 2 respectively
int counter=0;
vec betaG1null = betaG1;
betaG1null.subvec(start, end).fill(0);
vec betaG2null = betaG2;
betaG2null.subvec(start, end).fill(0);
vec rG1_null_fit = Xall_new.rows(0, Group1Visits_new*Q-1) * betaG1null;
vec rG2_null_fit = Xall_new.rows(Group1Visits_new*Q, Xall_new.n_rows-1) * betaG2null;
vec r_null_fit = join_cols(rG1_null_fit, rG2_null_fit);
mat r_null(Group1Visits_new + Group2Visits_new, Q, fill::zeros);
for(int v = 0; v < Group1Visits_new + Group2Visits_new; v++){
for(int q = 0; q < Q; q++){
r_null(v, q) = r_resid_new(v, q) + r_null_fit(counter);
counter ++;
}
}
List PsiBetaG1resid=getPsiBetaLongit(N1, Q, r_null.head_rows(Group1Visits_new),
RoyVars_new.head_slices(Group1Visits_new), step2_tol,
MaxIter, visvec_new.head(N1),
Xall_new.head_rows(Group1Visits_new*Q),
Psi0type, Psi1type, Group1Visits_new, Group1maxVis_new);
vec betaG1resid=PsiBetaG1resid("beta");
mat betaG1Varresid=PsiBetaG1resid("betaVar");
List PsiBetaG2resid=getPsiBetaLongit(N2, Q, r_null.tail_rows(Group2Visits_new),
RoyVars_new.tail_slices(Group2Visits_new), step2_tol,
MaxIter, visvec_new.tail(N2),
Xall_new.tail_rows(Group2Visits_new*Q),
Psi0type, Psi1type, Group2Visits_new, Group2maxVis_new);
vec betaG2resid=PsiBetaG2resid("beta");
mat betaG2Varresid=PsiBetaG2resid("betaVar");
double Tglobal = getStat(betaG1resid.subvec(start, end), betaG2resid.subvec(start, end),
betaG1Varresid.submat(start,start,end,end),
betaG2Varresid.submat(start,start,end,end));
int NPar=end-start+1;
vec Tlocal(NPar, fill::zeros);
for(int nPar = 0; nPar < NPar; nPar++){
Tlocal(nPar) = getStatUni(betaG1resid(start+nPar), betaG2resid(start+nPar),
betaG1Varresid(start+nPar, start+nPar),
betaG2Varresid(start+nPar, start+nPar));
}
return List::create(
_["Tglobal"] = Tglobal,
_["Tlocal"] = Tlocal
);
}
//function to permute data and calculate global main effect and global interaction test statistics for permuted data for longitudinal model
// [[Rcpp::export]]
List LongitTests(int N1, int N2, int Q, const mat& r_resid, const cube& RoyVars,
const mat& Xall, const vec& visvec, const vec& ids, double step2_tol,
string Psi0type, string Psi1type, int MaxIter,
const vec& betaG1, const vec& betaG2){
List permdat = Resample(N1+N2, Q, r_resid, RoyVars, Xall, visvec, false, ids);
cube RoyVars_new = permdat("RoyVars_new");
mat r_resid_new = permdat("r_new");
mat Xall_new = permdat("Xall_new");
vec visvec_new = permdat("visvec_new");
int Group1Visits_new = sum(visvec_new.head(N1));
int Group1maxVis_new = max(visvec_new.head(N1));
int Group2Visits_new = sum(visvec_new.tail(N2));
int Group2maxVis_new = max(visvec_new.tail(N2));
List Tmain = PermTestLongit(N1, N2, Q, r_resid_new, RoyVars_new, Xall_new, visvec_new,
step2_tol, Psi0type, Psi1type, MaxIter, betaG1,
betaG2, 0, Q-1, Group1Visits_new, Group2Visits_new,
Group1maxVis_new, Group2maxVis_new);
double TmainGlobal=Tmain("Tglobal");
vec TmainLocal=Tmain("Tlocal");
List Tint = PermTestLongit(N1, N2, Q, r_resid_new, RoyVars_new, Xall_new, visvec_new,
step2_tol, Psi0type, Psi1type, MaxIter, betaG1,
betaG2, Q, 2*Q-1, Group1Visits_new, Group2Visits_new,
Group1maxVis_new, Group2maxVis_new);
double TintGlobal=Tint("Tglobal");
vec TintLocal=Tint("Tlocal");
return List::create(
_["TmainGlobal"] = TmainGlobal,
_["TmainLocal"] = TmainLocal,