%% CodeForSWSStats.m % This file shows how to load the data contained in % volunteerSWSdata.txt into MATLAB, and relevant statistical tests % performed in the following manuscript: % Trutna Paley et al., "Repeatability of Rotational 3D Shear Wave % Elasticity Imaging Measurements in Skeletal Muscle" % For further information, contact Courtney Trutna Paley % courtney.trutna.paley@duke.edu % Code written in MATLAB R2022a %% Load data % load in data from CSV file to MATLAB table % Table should be 3600 x 7 A = readtable('dataset_subjectSWS.csv'); %% Seperate out data into MATLAB variables for easy processing % Plotting & statistical scripts have been implemented using matrix % variables with specific indexing. % Here the CSV file is loaded and converted to this format. % information for resizing arrays nvol=length(unique(A.SubjectNumber)); ntrials=length(unique(A.AcqTrialNum)); tmp=A.iSWSper3Dacq;tmp(isnan(tmp))=[]; % remove nans n2Dper3D=length(unique(tmp));clear tmp; % following should all be 10 subjects x 5 repeat grabs (size(X)=[10,5]) SWS_Along_3D_dom_singleloc=... reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,ntrials,nvol)'; SWS_Across_3D_dom_singleloc=... reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,ntrials,nvol)'; SWS_Along_3D_notdom_singleloc=... reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,ntrials,nvol)'; SWS_Across_3D_notdom_singleloc=... reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,ntrials,nvol)'; SWS_Along_3D_dom_reposition=... reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,ntrials,nvol)'; SWS_Across_3D_dom_reposition=... reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,ntrials,nvol)'; SWS_Along_3D_notdom_reposition=... reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,ntrials,nvol)'; SWS_Across_3D_notdom_reposition=... reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'3DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,ntrials,nvol)'; % Size of following should be 10 subjects x 5 repeat grabs x 8 2D-SWEI vals % per 3DSWEI acq (size(X)=[10,5,8]) SWSAlong_2D_dom_singleloc=permute(reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAcross_2D_dom_singleloc=permute(reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAlong_2D_notdom_singleloc=permute(reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAcross_2D_notdom_singleloc=permute(reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'singlelocation'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAlong_2D_dom_replace=permute(reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAcross_2D_dom_replace=permute(reshape(A(all([strcmp(A.LegType,'dominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAlong_2D_notdom_replace=permute(reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Along')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); SWSAcross_2D_notdom_replace=permute(reshape(A(all([strcmp(A.LegType,'nondominant'),... strcmp(A.AcqType,'repositioning'),... strcmp(A.SWEIType,'2DSWEI'),... strcmp(A.SWSDirection,'Across')],2),:).SWSvalue,n2Dper3D,ntrials,nvol),[3 2 1]); clear A nvol ntrials n2Dper3D %% Remove infinity values from 2D SWEI data % Shear wave speed estimation scripts sometimes output infinity for % 2D-SWEI values for certain error codes. These are error values, so % we will replace with NaN to avoid their inclusion in estimates of mean % or standard deviation SWSAlong_2D_dom_singleloc(isinf(SWSAlong_2D_dom_singleloc))=NaN; SWSAcross_2D_dom_singleloc(isinf(SWSAcross_2D_dom_singleloc))=NaN; SWSAlong_2D_notdom_singleloc(isinf(SWSAlong_2D_notdom_singleloc))=NaN; SWSAcross_2D_notdom_singleloc(isinf(SWSAcross_2D_notdom_singleloc))=NaN; SWSAlong_2D_dom_replace(isinf(SWSAlong_2D_dom_replace))=NaN; SWSAcross_2D_dom_replace(isinf(SWSAcross_2D_dom_replace))=NaN; SWSAlong_2D_notdom_replace(isinf(SWSAlong_2D_notdom_replace))=NaN; SWSAcross_2D_notdom_replace(isinf(SWSAcross_2D_notdom_replace))=NaN; %% Calculate coefficent of variation for each measurement, each leg % intermediate values used in next section CV_Along_dom_singleloc=std(SWS_Along_3D_dom_singleloc,[],2)./mean(SWS_Along_3D_dom_singleloc,2); CV_Across_dom_singleloc=std(SWS_Across_3D_dom_singleloc,[],2)./mean(SWS_Across_3D_dom_singleloc,2); CV_Along_notdom_singleloc=std(SWS_Along_3D_notdom_singleloc,[],2)./mean(SWS_Along_3D_notdom_singleloc,2); CV_Across_notdom_singleloc=std(SWS_Across_3D_notdom_singleloc,[],2)./mean(SWS_Across_3D_notdom_singleloc,2); CV_Along_dom_replace=std(SWS_Along_3D_dom_reposition,[],2)./mean(SWS_Along_3D_dom_reposition,2); CV_Across_dom_replace=std(SWS_Across_3D_dom_reposition,[],2)./mean(SWS_Across_3D_dom_reposition,2); CV_Along_notdom_replace=std(SWS_Along_3D_notdom_reposition,[],2)./mean(SWS_Along_3D_notdom_reposition,2); CV_Across_notdom_replace=std(SWS_Across_3D_notdom_reposition,[],2)./mean(SWS_Across_3D_notdom_reposition,2); CV2D_Along_dom_replace=std(SWSAlong_2D_dom_replace,[],[2 3],'omitnan')./mean(SWSAlong_2D_dom_replace,[2 3],'omitnan'); CV2D_Across_dom_replace=std(SWSAcross_2D_dom_replace,[],[2 3],'omitnan')./mean(SWSAcross_2D_dom_replace,[2 3],'omitnan'); CV2D_Along_notdom_replace=std(SWSAlong_2D_notdom_replace,[],[2 3],'omitnan')./mean(SWSAlong_2D_notdom_replace,[2 3],'omitnan'); CV2D_Across_notdom_replace=std(SWSAcross_2D_notdom_replace,[],[2 3],'omitnan')./mean(SWSAcross_2D_notdom_replace,[2 3],'omitnan'); %% Calculate average CVs for directions along and across, different acquisition types % These are the values reported in Table II averageCV_singlelocation_along=mean([CV_Along_dom_singleloc(:); CV_Along_notdom_singleloc(:)]) averageCV_singlelocation_across=mean([CV_Across_dom_singleloc(:); CV_Across_notdom_singleloc(:)]) averageCV_repositioning3D_along=mean([CV_Along_dom_replace(:); CV_Along_notdom_replace(:)]) averageCV_repositioning3D_across=mean([CV_Across_dom_replace(:); CV_Across_notdom_replace(:)]) averageCV_repositioning2D_along=mean([CV2D_Along_dom_replace(:); CV2D_Along_notdom_replace(:)]) averageCV_repositioning2D_across=mean([CV2D_Across_dom_replace(:); CV2D_Across_notdom_replace(:)]) %% Calculate Example Confidence Intervals % calculate standard deviations averagestd_repositioning_along=mean([std(SWS_Along_3D_dom_reposition,[],2); std(SWS_Along_3D_notdom_reposition,[],2)]); averagestd_repositioning_across=mean([std(SWS_Across_3D_dom_reposition,[],2); std(SWS_Across_3D_notdom_reposition,[],2)]); averagestd_2Drepositioning_along=mean([std(SWSAlong_2D_dom_replace,[],[2,3],'omitnan'); std(SWSAlong_2D_notdom_replace,[],[2,3])],'omitnan'); averagestd_2Drepositioning_across=mean([std(SWSAcross_2D_dom_replace,[],[2,3],'omitnan'); std(SWSAcross_2D_notdom_replace,[],[2,3])],'omitnan'); % Assuming our standard deviation (s) is the same as what would be seen in % a new hypothetical study with n measurements, % the margin of error is Z*s/sqrt(n). Z=1.96; %95% confidence nvectortest=[1 5 10 25]; % test these options of n %Confidence intervals: values reported in Table III. alongCIs_3D=Z*averagestd_repositioning_along./sqrt(nvectortest) acrossCIs_3D=Z*averagestd_repositioning_across./sqrt(nvectortest) alongCIs_2D=Z*averagestd_2Drepositioning_along./sqrt(nvectortest) acrossCIs_2D=Z*averagestd_2Drepositioning_across./sqrt(nvectortest) %% Perform Leg Dominance Hypothesis Testing % test normality of differences for paired ttest % shapiro wilks test coded by Ahmed Ben Saïda, available: % https://www.mathworks.com/matlabcentral/fileexchange/13964-shapiro-wilk-and-shapiro-francia-normality-tests try [H_NormalAlong, pValue_NormalAlong, ~] = swtest(mean(SWS_Along_3D_dom_reposition,2)-mean(SWS_Along_3D_notdom_reposition,2)) [H_NormalAcross, pValue_NormalAcross, ~] = swtest(mean(SWS_Across_3D_dom_reposition,2)-mean(SWS_Across_3D_notdom_reposition,2)) catch warning('Cannot complete normality testing. Please add swtest() to MATLAB path') end % Test along and across the fibers seperately using paired t test. % please note the ttest used below requires the % Statistics and Machine Learning Toolbox [H_DifferentAlong,p_DifferentAlong]=ttest(mean(SWS_Along_3D_dom_reposition,2),mean(SWS_Along_3D_notdom_reposition,2)) % paired T-test on dom vs not dom [H_DifferentAcross,p_DifferentAcross]=ttest(mean(SWS_Across_3D_dom_reposition,2),mean(SWS_Across_3D_notdom_reposition,2)) % paired T-test on dom vs not dom