version 0.1.0

DESCRIPTION

@@ -0,0 +1,29 @@
+Package: dfConn
+Type: Package
+Title: Dynamic Functional Connectivity Analysis
+Version: 0.1.0
+Authors@R: c(person("Zikai", "Lin", email = "ziklin@iu.edu", role = c("aut", "cre")),
+            person("Maria", "Kudela", email = "maria.kudela@gmail.com", role = c("aut")),
+            person("Jaroslaw","Harezlak",email = "harezlak@iu.edu", role = c("aut")),
+            person("Mario","Dzemidzic", email = "mdzemidz@iupui.edu", role = c("aut")))
+Maintainer: Zikai Lin <ziklin@iu.edu>
+Description: An implementation of multivariate linear process bootstrap (MLPB) method and sliding window technique to assess the dynamic functional connectivity (dFC) estimate by providing its confidence bands, based on Maria Kudela (2017) <doi: 10.1016/j.neuroimage.2017.01.056>.
+  It also integrates features to visualize non-zero coverage for selected a-priori regions of interest estimated by the dynamic functional connectivity model (dFCM) and dynamic functional connectivity (dFC) curves for reward-related a-priori regions of interest where the activation-based analysis reported.
+License: MIT + file LICENSE
+Encoding: UTF-8
+LazyData: true
+Depends: R (>= 2.10)
+Suggests: iterators, testthat, itertools
+Imports: doParallel, nlme, parallel, foreach, ggplot2, knitr,
+        RColorBrewer, fields, latex2exp, mgcv, gplots, splines, stats,
+        stringr, graphics, data.table, gtools, Rcpp (>= 0.12.18)
+LinkingTo: Rcpp, RcppArmadillo
+RoxygenNote: 6.1.1
+NeedsCompilation: yes
+Packaged: 2019-03-22 16:22:03 UTC; ziklin
+Author: Zikai Lin [aut, cre],
+  Maria Kudela [aut],
+  Jaroslaw Harezlak [aut],
+  Mario Dzemidzic [aut]
+Repository: CRAN
+Date/Publication: 2019-03-23 09:13:24 UTC

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2019
+COPYRIGHT HOLDER: Lin, Zikai

MD5

@@ -0,0 +1,40 @@
+fa28a1cdc5b5805cbfec838ac043fc4d *DESCRIPTION
+8dc4dd367716da3d741838c6a8e830ed *LICENSE
+93b7a7047065da580135d8dce6890273 *NAMESPACE
+1e70ac2c81d1c915a6bf6206f4b5638c *R/MLPB3.R
+09467175ccefa10462de4608b2b83769 *R/RcppExports.R
+c892bfcedfa1dc4dc17d51670001314d *R/data_summary_coverage.r
+771904ec595c41c37a5eaaa284955008 *R/dfConn.R
+66199f7cbe55168ee969ec3faeb4c8ac *R/lme_model0.r
+280e7ef9a179406017ccad79de2fc7ac *R/lme_model6.r
+1ad2e4f3571ca0d04394a2ca5ee582f6 *R/lmm_ci_vis.R
+e7420249fc19c4e3cb6791a217683ac1 *R/mlpb.R
+48345987ba7c60695901db2d4eb808ff *R/nzc_mat_construct.r
+a2ca620d8352fc8a38d8c2ebbf8e4d48 *R/nzc_mat_vis.r
+a2f439f6db7c754a7f0182f6143947f4 *R/plots.r
+61f720d5cb54596384caeae466f95993 *R/setters.R
+8d39b9c394ee48e22b7e392e9db4d594 *R/summary_table.r
+852948da1bed1d808d2a3739de5ae977 *R/sysdata.rda
+1137d1ad613ebc8fba702abfb16ca33d *R/utils.r
+929fb0a54af8e643467cf9a4a5a63a25 *README.md
+1dc116c571fd89cef1cbde6dadb83f66 *data/DynModel_results.rda
+c48f1474670b3f46d96982ecdade3d4d *data/MLPB_output_median.rda
+24944c332629f20dd1c739319178572f *data/coverage.tbl.list.rda
+d4fbd1314ed2648a8021c1336de9e6c1 *data/fMRI_dataList.rda
+c9bb9b3e1e798019ffa7f7483640cd0e *data/fMRI_dataList_shrinked.rda
+5073b53fa8f3088d8a6e4803a1f5e692 *man/DynModel_results.Rd
+709765a2b7c36d468f0676dfe24116f6 *man/MLPB3.Rd
+52d620f6e47149e5ab257c5ab674e824 *man/MLPB_boot.Rd
+ab420afbcb3b38f9209756efe02161a7 *man/MLPB_output_median.Rd
+66fb1817a2bcbca31c12fcb78a571dda *man/coverage.tbl.list.Rd
+625cd173de4cf8b896e448d235bb5225 *man/data_summary_coverage.Rd
+0a987d54dde8932b4799ed5867749d36 *man/fMRI_dataList.Rd
+97f3abf819c842e4d85af868b90cbd20 *man/fMRI_dataList_shrinked.Rd
+c1d5ee13d0e7343a822a3a1cdd63cb9d *man/lmmConn.Rd
+f8041479fe05ef1ed078a5af8654a54a *man/lmmDyn.Rd
+294039488892ec38f4342b58d921a1ca *man/nzc_vis.Rd
+f103a402addfdebb68cca4c21b339339 *man/plot_generate_timeDependence.Rd
+23b3dd69cab930e98293c85d89100dfe *src/MLPB3.cpp
+61059660eb073d93e00e8ee054237071 *src/Makevars
+61059660eb073d93e00e8ee054237071 *src/Makevars.win
+af378eca0bc4b9114d08e0c7d121036b *src/RcppExports.cpp

NAMESPACE

@@ -0,0 +1,21 @@
+# Generated by roxygen2: do not edit by hand
+
+export(MLPB3)
+export(MLPB_boot)
+export(data_summary_coverage)
+export(lmmConn)
+export(lmmDyn)
+export(nzc_vis)
+export(plot_generate_timeDependence)
+import(Rcpp)
+importFrom(Rcpp,evalCpp)
+importFrom(data.table,rbindlist)
+importFrom(gtools,combinations)
+importFrom(parallel,makePSOCKcluster)
+importFrom(stats,cor)
+importFrom(stats,mad)
+importFrom(stats,median)
+importFrom(stats,qnorm)
+importFrom(stats,quantile)
+importFrom(utils,write.csv)
+useDynLib(dfConn)

R/MLPB3.R

@@ -0,0 +1,230 @@
+#' @title Multivariate Linear Bootstrapping - core function
+#' @param X data matrix
+#' @param Boot integer, number of bootstrap samples to be generated
+#' @param l numeric, the banding parameter, default is 1
+#' @param eps numeric the parameters for making Gamma_kappa_matrix positive definite if necessary, default is 1;
+#' @param beta numeric the parameters for making Gamma_kappa_matrix positive definite if necessary, default is 1;
+#' @param l_automatic numeric the banding parameter, default is 1, data-adaptively;
+#' @param l_automatic_local numeric, the banding parameter default is 0 
+#' @return  a numeric matrix with n*boot rows and m columns, where n, m refer to the number of rows and number of columns in input data matrix, and boot refer to number of bootstrap samples to be generated.
+#' @examples 
+#' # Multivariate linear bootstrapping on a random matrix with 2 rows and 4 columns
+#' MLPB3(matrix(rnorm(16),2,4), 3)
+#' 
+#' @export
+#' @useDynLib dfConn
+
+
+MLPB3 <- function(X, Boot, l = 1, eps = 1, beta = 1, l_automatic = 0, l_automatic_local = 0) {
+    requireNamespace("stats")
+    # Dimension of observed values
+    n1 <- dim(X)[1]
+
+    # sample sizes of (n1-dimensional) data observed
+    n2 <- dim(X)[2]
+    N <- n1 * n2
+    # Boot<-100
+    l_matrix <- matrix(rep(l, n1^2), n1, n1)
+
+    # Defining function kappa (trapezoid [cf. McMurry & Politis (2010)])
+    kappa <- function(x, l) {
+        kappa <- matrix(0, dim(l)[1], dim(l)[2])
+        kappa[l > 0] <- 2 - abs(x/l[l > 0])
+        kappa[abs(x/l) <= 1 & l != 0] <- 1
+        kappa[abs(x/l) > 2 & l != 0] <- 0
+        kappa
+    }
+
+    # Rearranging the data columnwise in one large n1*n2-dimensional vector
+    X_vec <- X
+    dim(X_vec) <- NULL
+
+    ######################################################################## Here starts the computation of covariances ###
+
+    # Estimating all multivariate covariances and putting them in one large n1*(2n2-1) matrix C_matrix<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1)
+    # C_matrix_kappa_trans<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1) R_matrix_kappa_trans<-matrix(rep(0,(n1*n1*(2*n2-1))),n1,(2*n2-1)*n1)
+    C_matrix <- C_matrix_kappa_trans <- R_matrix_kappa_trans <- matrix(0, n1, (2 * n2 - 1) * n1)
+
+    # Sample mean
+    if (n1 == 1) 
+        X_mean <- mean(X) else X_mean <- rowMeans(X)
+    # if(n1!=1){ X_mean<-rep(0,n1) for(k in 1:n1){ X_mean[k]<-mean(X[k,]) } }
+
+    # Computing C_matrix
+    for (k in (-(n2 - 1)):(n2 - 1)) {
+        blub <- 0
+        if (n1 == 1) 
+            C_matrix[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- (1/n2) * t(X[(max(1, 1 - k) + k):(min(n2, n2 - k) + k)] - X_mean) %*% (X[(max(1, 1 - k)):(min(n2, 
+                n2 - k))] - X_mean) else C_matrix[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- (1/n2) * (X[, (max(1, 1 - k) + k):(min(n2, n2 - k) + k)] - X_mean) %*% t(X[, (max(1, 1 - 
+            k)):(min(n2, n2 - k))] - X_mean)
+        # if(n1!=1) C_matrix[,((((n2+k)-1)*n1+1):(((n2+k)-1)*n1+n1))]<-(1/n2)*(X[,(max(1,1-k)+k):(min(n2,n2-k)+k)]-X_mean)%*%t(X[,(max(1,1-k)):(min(n2,n2-k))]-X_mean)
+    }
+
+    if (l_automatic == 1) {
+
+        ################################## Computing Autocorrelations ###
+
+        R_matrix <- matrix(0, n1, (2 * n2 - 1) * n1)
+        denominator_matrix <- matrix(0, n1, n1)
+        diag(denominator_matrix) <- sqrt(diag(C_matrix[, (n1 * (n2 - 1) + 1):(n1 * n2)]))
+        denominator_matrix <- solve(denominator_matrix)
+        for (k in 1:(2 * n2 - 1)) {
+            R_matrix[, ((k - 1) * n1 + 1):(k * n1)] <- denominator_matrix %*% C_matrix[, ((k - 1) * n1 + 1):(k * n1)] %*% denominator_matrix
+        }
+
+        M_zero <- 2
+        M_zero_sqrt <- M_zero * sqrt(log10(n2)/n2)
+
+        if (l_automatic_local == 0) {
+            h <- 1
+            marker <- 0
+            while (marker == 0) {
+                # Kn_vec<-seq(1,max(5,log10(n2))) if(max((abs(R_matrix[,((2*n2-1)+h*n1):((2*n2-1+n1-1)+n1*(max(5,log10(n2))+h-1))])))<M_zero*sqrt(log10(n2)/n2)) marker<-1
+                if (max((abs(R_matrix[, (n1 * (n2 - 1) + 1 + h * n1):(n1 * n2 + n1 * (max(5, log10(n2)) + h - 1))]))) < M_zero_sqrt) 
+                  marker <- 1
+                l_adapt <- h - 1
+                h <- h + 1
+            }
+            l <- l_adapt
+            l_matrix_adapt <- matrix(rep(l, n1^2), n1, n1)
+            l_matrix <- l_matrix_adapt
+        }
+
+        if (l_automatic_local == 1) {
+            l_matrix_adapt <- matrix(0, n1, n1)
+            marker <- matrix(0, n1, n1)
+            for (j in 1:n1) {
+                for (k in 1:n1) {
+                  h <- 1
+                  while (marker[j, k] == 0) {
+                    if (max(abs(R_matrix[j, seq((n1 * (n2 - 1) + 1 + h * n1 + k - 1), (n1 * n2 + n1 * (max(5, log10(n2)) + h - 1) + k - 1), by = n1)])) < M_zero_sqrt) 
+                      marker[j, k] <- 1
+                    l_matrix_adapt[j, k] <- h - 1
+                    h <- h + 1
+                  }
+                }
+            }
+            l_matrix <- l_matrix_adapt
+        }
+
+    }
+
+    # Computing C_matrix_kappa_trans
+    for (k in (-(n2 - 1)):(n2 - 1)) {
+        blub <- 0
+        if (n1 == 1) 
+            C_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- (kappa(k, l_matrix)/n2) * t(X[(max(1, 1 - k) + k):(min(n2, n2 - k) + k)] - X_mean) %*% 
+                (X[(max(1, 1 - k)):(min(n2, n2 - k))] - X_mean)
+        if (n1 != 1 && k < 0) 
+            C_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- t(t((kappa(k, l_matrix)/n2)) * (X[, (max(1, 1 - k) + k):(min(n2, n2 - k) + k)] - 
+                X_mean) %*% t(X[, (max(1, 1 - k)):(min(n2, n2 - k))] - X_mean))
+        if (n1 != 1 && k == 0) 
+            C_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- t((1/n2) * (X[, (max(1, 1 - k) + k):(min(n2, n2 - k) + k)] - X_mean) %*% t(X[, 
+                (max(1, 1 - k)):(min(n2, n2 - k))] - X_mean))
+        if (n1 != 1 && k > 0) 
+            C_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- t((kappa(k, l_matrix)/n2) * (X[, (max(1, 1 - k) + k):(min(n2, n2 - k) + k)] - 
+                X_mean) %*% t(X[, (max(1, 1 - k)):(min(n2, n2 - k))] - X_mean))
+    }
+    ############################################################################# Computing the diagonal matrix for standardization to get autocorrelations (as above)
+    denominator_matrix <- matrix(0, n1, n1)
+    diag(denominator_matrix) <- sqrt(diag(C_matrix[, (n1 * (n2 - 1) + 1):(n1 * n2)]))
+    denominator_matrix <- solve(denominator_matrix)
+
+    # Computing R_matrix_kappa_trans (for equivariance to construct positive definite estimator)
+    for (k in (-(n2 - 1)):(n2 - 1)) {
+        blub <- 0
+        if (n1 == 1) 
+            R_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- (1/stats::var(as.vector(X))) * (kappa(k, l_matrix)/n2) * t(X[(max(1, 1 - k) + 
+                k):(min(n2, n2 - k) + k)] - X_mean) %*% (X[(max(1, 1 - k)):(min(n2, n2 - k))] - X_mean)
+        if (n1 != 1 && k < 0) 
+            R_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- denominator_matrix %*% t(t((kappa(k, l_matrix)/n2)) * (X[, (max(1, 1 - k) + 
+                k):(min(n2, n2 - k) + k)] - X_mean) %*% t(X[, (max(1, 1 - k)):(min(n2, n2 - k))] - X_mean)) %*% denominator_matrix
+        if (n1 != 1 && k == 0) 
+            R_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- denominator_matrix %*% t((1/n2) * (X[, (max(1, 1 - k) + k):(min(n2, n2 - k) + 
+                k)] - X_mean) %*% t(X[, (max(1, 1 - k)):(min(n2, n2 - k))] - X_mean)) %*% denominator_matrix
+        if (n1 != 1 && k > 0) 
+            R_matrix_kappa_trans[, ((((n2 + k) - 1) * n1 + 1):(((n2 + k) - 1) * n1 + n1))] <- denominator_matrix %*% t((kappa(k, l_matrix)/n2) * (X[, (max(1, 1 - k) + k):(min(n2, 
+                n2 - k) + k)] - X_mean) %*% t(X[, (max(1, 1 - k)):(min(n2, n2 - k))] - X_mean)) %*% denominator_matrix
+    }
+
+    # Arranging the autocovariance matrices in R_matrix_kappa_trans in one large matrix corresponding to X_vec
+    R_kappa_matrix <- matrix(rep(0, N * N), N, N)
+    for (i in 1:n2) {
+        R_kappa_matrix[(((i - 1) * n1 + 1):((i - 1) * n1 + n1)), ] <- R_matrix_kappa_trans[, (n1 * (n2 - i) + 1):(n1 * (2 * n2 - i))]
+    }
+
+    R_kappa_matrix_hilf <- R_kappa_matrix
+
+    # Computing the eigenvalues and eigenvectors of Gamma_kappa_matrix
+    spectraldecomposition <- eigen(R_kappa_matrix)
+
+    # orthogonal matrix
+    S <- spectraldecomposition$vectors
+
+    # diagonal matrix
+    D <- diag(spectraldecomposition$values)
+
+    D_eps_counter <- 0
+    # Manipulating the diagonal matrix D to ensure positive definiteness
+    if (min(diag(D)) <= eps * N^(-beta)) {
+        D_eps_counter <- D_eps_counter + 1
+        D_eps <- D
+        D_non_definite <- D
+        # for(k in 1:N){ D_eps[k,k]<-max(D[k,k],eps*N^(-beta)) }
+        diag(D_eps) <- ifelse(diag(D) > eps * N^(-beta), diag(D), eps * N^(-beta))
+        # Computing the banded positive definite covariance matrix estimator
+        R_kappa_matrix <- S %*% D_eps %*% t(S)
+    }
+
+    # Get back to autocovariances!
+    Gamma_kappa_matrix <- matrix(rep(0, N * N), N, N)
+    nominator_matrix <- solve(denominator_matrix)
+    for (i in 1:n2) {
+        for (j in 1:n2) {
+            Gamma_kappa_matrix[(((i - 1) * n1 + 1):((i - 1) * n1 + n1)), (((j - 1) * n1 + 1):((j - 1) * n1 + n1))] <- nominator_matrix %*% R_kappa_matrix[(((i - 1) * n1 + 
+                1):((i - 1) * n1 + n1)), (((j - 1) * n1 + 1):((j - 1) * n1 + n1))] %*% nominator_matrix
+        }
+    }
+
+
+    # Computing the (lower-left) Cholesky decomposition
+    L <- t(chol(Gamma_kappa_matrix))
+    L_inv <- solve(L)
+
+    # Multilpication of L and centered X_vec that is X_vec_cent X_cent<-matrix(rep(0,n1*n2),n1,n2) for(k in 1:n1){ X_cent[k,]<-X[k,]-X_mean[k] }
+    X_cent <- X - rowMeans(X)
+
+    # Rearranging the data columnwise in one large n1*n2-dimensional vector
+    X_vec_cent <- X_cent
+    dim(X_vec_cent) <- NULL
+
+    # 'Whitening' the centered data
+    W_vec <- L_inv %*% X_vec_cent
+
+    # Computing the centered and standardized residuals Z_vec
+    W_vec_cent <- W_vec - mean(W_vec)
+    W_vec_var <- stats::var(W_vec)
+    dim(W_vec_var) <- NULL
+    Z_vec <- (1/sqrt(W_vec_var)) * W_vec_cent
+
+    W_matrix <- matrix(W_vec, n1, n2)
+    # W_matrix_cent<-matrix(rep(0,n1*n2),n1,n2) for(k in 1:n1){ W_matrix_cent[k,]<-W_matrix[k,]-mean(W_matrix[k,]) }
+
+    W_matrix_cent <- W_matrix - rowMeans(W_matrix)
+
+    W_matrix_var <- (1/n2) * W_matrix_cent %*% t(W_matrix_cent)
+    W_matrix_var_inv <- solve(t(chol(W_matrix_var)))
+
+    # W_matrix_stud<-matrix(rep(0,n1*n2),n1,n2) for(k in 1:n2){ W_matrix_stud[,k]<-W_matrix_var_inv%*%W_matrix_cent[,k] }
+    W_matrix_stud <- W_matrix_var_inv %*% W_matrix_cent
+
+    Z_vec <- W_matrix_stud
+    dim(Z_vec) <- NULL
+
+    ########################################################################## Here starts the LPB bootstrap loop ###
+
+
+    # rewrite LPB bootstrap loop using Rcpp, Xiaochun Li
+    sv_draw <- replicate(Boot, sample(1:n2, n2, replace = TRUE))
+    boot_func(Boot, sv_draw, matrix(Z_vec, nrow = n1), L)
+}

R/RcppExports.R

@@ -0,0 +1,27 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+myseq <- function(start, end, by) {
+    .Call('_dfConn_myseq', PACKAGE = 'dfConn', start, end, by)
+}
+
+boot_func <- function(Boot, sw_mat, Z_mat, L) {
+    .Call('_dfConn_boot_func', PACKAGE = 'dfConn', Boot, sw_mat, Z_mat, L)
+}
+
+mm_func <- function(Xjj, F, mm, n_boot, window_size = 20L, n_sig = 750L) {
+    .Call('_dfConn_mm_func', PACKAGE = 'dfConn', Xjj, F, mm, n_boot, window_size, n_sig)
+}
+
+boot_cjj_func <- function(Xjj, MLPB3, boot_rep = 250L, n_boot = 1L, n_sig = 750L, window_size = 20L) {
+    .Call('_dfConn_boot_cjj_func', PACKAGE = 'dfConn', Xjj, MLPB3, boot_rep, n_boot, n_sig, window_size)
+}
+
+a_func <- function(bootjj, up_limit = 731L, n_sig = 750L, window_size = 20L) {
+    .Call('_dfConn_a_func', PACKAGE = 'dfConn', bootjj, up_limit, n_sig, window_size)
+}
+
+bootcjj <- function(X, MLPB3, up_limit = 731L, window_size = 20L, boot_rep = 250L, n_boot = 1L, n_sig = 750L) {
+    .Call('_dfConn_bootcjj', PACKAGE = 'dfConn', X, MLPB3, up_limit, window_size, boot_rep, n_boot, n_sig)
+}
+