%========================================================================== % NMR DATA - VEHICLE PAPER %========================================================================== % Load data load CINMR_RAW.mat; %-------------------------------------------------------------------------- % DATA PREPARATION %-------------------------------------------------------------------------- % Set detection limit Raw1=Raw; Raw1(Raw<2000000)=0; % Normalize wrt Creatinine C1=find(strcmp(Bins,'B_4.07')==1); C2=find(strcmp(Bins,'B_4.05')==1); C3=find(strcmp(Bins,'B_3.05')==1); CREAT=Raw(:,C1) + Raw(:,C2); Data=Raw1./repmat(CREAT,1,size(Raw1,2)); Data1=Data; Data1(:,[C1, C2, C3])=[]; Bins1=Bins; Bins1([C1, C2, C3])=[]; %-------------------------------------------------------------------------- % QC ANALYSIS %-------------------------------------------------------------------------- iQC=find(strcmp(Treatment,'QC')==1); QCData=center(log(Data1(iQC,:)+1)); % Perform PCA using PLS Toolbox PCAOutliers(QCData,'QC',0.9,num2cell(Batch(iQC)),3,''); % Evaluate scores and Hotteling's T plots to identify & remove outlying batches idrop=find(or(Batch==2,Batch==3)==1); Data2=Data1; Data2(idrop,:)=[]; TR=Treatment; TR(idrop)=[]; T=Time; T(idrop)=[]; B=Batch; B(idrop)=[]; %-------------------------------------------------------------------------- % DATA PREPARATION %-------------------------------------------------------------------------- % Select vehicle data i2use=find(strcmp(TR,'Vehicle')==1); USE=Data2(i2use,:); TR=TR(i2use); T=T(i2use); B=B(i2use); % Remove Variables with 50%+ zeroes... dat_1=USE(T==-1,:); dat0= USE(T==0,:); dat1= USE(T==1,:); dat2= USE(T==2,:); dat3= USE(T==3,:); dat4= USE(T==4,:); USE2 = []; Bins2 = {}; outmets={}; z=0.5; for c = 1:size(USE,2); if ((sum(dat_1(:,c)==0)/size(dat_1,1))>z)&&... ((sum(dat0(:,c)==0)/size(dat0,1))>z)&&... ((sum(dat1(:,c)==0)/size(dat1,1))>z)&&... ((sum(dat2(:,c)==0)/size(dat2,1))>z)&&... ((sum(dat3(:,c)==0)/size(dat3,1))>z)&&... ((sum(dat4(:,c)==0)/size(dat4,1))>z); outmets = [outmets; Bins1(c)]; else USE2 = [USE2 , USE(:,c)]; Bins2 = [Bins2; Bins1(c)]; end; end; % Scale USE3=center(log(USE2+1)); % Reset Batch numbers fixB=B; fixB(and(B>1,B<8)==1)=B(and(B>1,B<8)==1)-2; fixB(B>8)=B(B>8)-3; %-------------------------------------------------------------------------- % CROSS-SECTIONAL %-------------------------------------------------------------------------- % The Compare2 function calculates: % i) Univariate statistics using the statistics toolbox % ii) PCA using the PLS Toolbox % iii) PLS-DA using the PLS Toolbox i2use=find(or(T==-1,T==0)==1); t1='-1'; t2=' 0'; Compare2(USE2(i2use,:),USE3(i2use,:),USE2(i2use,:),... str2cell(num2str(T(i2use))),... t1,t2,'b','c','s','o',... str2cell(num2str(B(i2use))),... Bins2','Vehicle',... 0.9,0,'yes','between',1,1); i2use=find(or(T==0,T==1)==1); t1='0'; t2='1'; Compare2(USE2(i2use,:),USE3(i2use,:),USE2(i2use,:),... str2cell(num2str(T(i2use))),... t1,t2,'b','c','s','o',... str2cell(num2str(B(i2use))),... Bins2','Vehicle',... 0.9,0,'yes','between',1,1); i2use=find(or(T==0,T==2)==1); t1='0'; t2='2'; Compare2(USE2(i2use,:),USE3(i2use,:),USE2(i2use,:),... str2cell(num2str(T(i2use))),... t1,t2,'b','c','s','o',... str2cell(num2str(B(i2use))),... Bins2','Vehicle',... 0.9,0,'yes','between',1,1); i2use=find(or(T==0,T==3)==1); t1='0'; t2='3'; Compare2(USE2(i2use,:),USE3(i2use,:),USE2(i2use,:),... str2cell(num2str(T(i2use))),... t1,t2,'b','c','s','o',... str2cell(num2str(B(i2use))),... Bins2','Vehicle',... 0.9,0,'yes','between',1,1); i2use=find(or(T==0,T==4)==1); t1='0'; t2='4'; Compare2(USE2(i2use,:),USE3(i2use,:),USE2(i2use,:),... str2cell(num2str(T(i2use))),... t1,t2,'b','c','s','o',... str2cell(num2str(B(i2use))),... Bins2','Vehicle',... 0.9,0,'yes','between',1,1); %-------------------------------------------------------------------------- % Repeated measures ANOVA %-------------------------------------------------------------------------- % The RMAOV1_mvr function performs one-way repeated measures ANOVA and % is a slight variation of the original RMAOV1 function: % Trujillo-Ortiz, A., R. Hernandez-Walls and R.A. Trujillo-Perez. (2004). % RMAOV1:One-way repeated measures ANOVA. % A MATLAB file. [WWW document]. % URL http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=5576 DS=[{'Repeated Measures ANOVA p-values'},{'_'},{'_'}]; DS=[DS;[{'Variable'},{'Subject Effect'},{'Time Effect'}]]; Psubject=[]; Ptreatment=[]; Ptime=[]; Pinteraction=[]; for i=1:length(Bins2); X=[USE2(:,i),T+2,fixB]; [pS,pT] = RMAOV1_mvr(X,0.05); DS=[DS;[Bins2(i),num2cell(pS),num2cell(pT)]]; Psubject=[Psubject;pS]; Ptime=[Ptime;pT]; end xlswrite('Repeated Measures ANOVA Vehicle.xlsx', DS, 'P-values', 'A1'); %-------------------------------------------------------------------------- % ASCA %-------------------------------------------------------------------------- % The asca function performs principal component analysis on the % level averages of each factor, as provided by: % Gooitzen Zwanenburg (Copyright) % This code is available under APACHE Licence 2.0 % http://www.apache.org/licenses/LICENSE-2.0.html ASCA_V = asca(USE3, [fixB,T+2],[], 1); close all; % ASCA "VIPs" % Sum of squared loadings loads=ASCA_V.factors.loadings{2,1}; % TIME FACTOR SSQL=sum((loads(:,1:2).^2)'); eigvals=ASCA_V.factors.singular{2,1}; VarExpl = eigvals./sum(eigvals); DS=[{'Variable'},{'Sum of Squared Loadings for First 2 LVs'}]; DS=[DS;[Bins2,num2cell(SSQL)']]; xlswrite('VIP ASCA TIME FACTOR FOR VEHICLE.xlsx', DS, 'Variable SSQL List', 'A1'); %-------------------------------------------------------------------------- % PCA UNFOLDED %-------------------------------------------------------------------------- % UNFOLDING - NB MUST BE ORDERED BY SUBJECT unfold=[]; unfold=[unfold,USE3(T==-1,:)]; unfold=[unfold,USE3(T==0,:)]; unfold=[unfold,USE3(T==1,:)]; unfold=[unfold,USE3(T==2,:)]; unfold=[unfold,USE3(T==3,:)]; unfold=[unfold,USE3(T==4,:)]; % PCA % Perform PCA using the PLS Toolbox options = pca('options'); options.preprocessing = {[] []}; modelpca = pca(unfold,2); loadings = modelpca.loads{2,1}; scores = modelpca.loads{1,1}; eigvals = modelpca.detail.ssq(:,3); VarExpl = eigvals./sum(eigvals); close all;