too many arguments error using uncertainty and sensitivity code

Please, I need help on how to carry out uncertainty and sensitivity analysis of a disease model. I came acorss PRCC code in this link. http://malthus.micro.med.umich.edu/lab/usadata

I try to run the code on my laptop, It showed many error arguments. Also, I do not know which one of the codes run first.

For the code " ODE_LHS "

function dydt=ODE_LHS(t,y,LHSmatrix,x,runs)
%% PARAMETERS %%
Parameter_settings_LHS;

s=LHSmatrix(x,1);
muT=LHSmatrix(x,2);
r=LHSmatrix(x,3);
k1=LHSmatrix(x,4);
k2=LHSmatrix(x,5);
mub=LHSmatrix(x,6);
N=LHSmatrix(x,7);
muV=LHSmatrix(x,8);
dummy=LHSmatrix(x,9);

% [T] CD4+ uninfected: Tsource + Tprolif - Tinf
Tsource = s - muT*y(1);
Tprolif = r*y(1)*(1-(y(1)+y(2)+y(3))/Tmax);
Tinf = k1*y(1)*y(4);

% [T1] CD4+ latently infected: Tinf - T1death - T1inf
T1death = muT*y(2);
T1inf = k2*y(2);

% [T2] CD4+ actively infected: T1inf - T2death
T2death = mub*y(3);

% [V] Free infectious virus: Vrelease - Tinf - Vdeath
Vrelease = N*T2death;
Vdeath = muV*y(4);

dydt = [Tsource + Tprolif - Tinf;
        Tinf - T1death - T1inf;
        T1inf - T2death;
        Vrelease - Tinf - Vdeath];

I got this error message

??? Input argument "x" is undefined.

Error in ==> ODE_LHS at 5

s=LHSmatrix(x,1);

Also, for this follwing code;

function s=LHS_Call(xmin,xmean,xmax,xsd,nsample,distrib,threshold)
% s=latin_hs(xmean,xsd,nsample,nvar)
% LHS from normal distribution, no correlation
% method of Stein
% Stein, M. 1987. Large Sample Properties of Simulations Using Latin Hypercube Sampling. 
%                 Technometrics 29:143-151

if nsample==1
    s=xmean;
    return
end
if nargin<7
    threshold=1e20;
end

[sample,nvar]=size(xmean);
if distrib == 'norm'   % you only need to specify xmean & xsd
    ran=rand(nsample,nvar);
    s=zeros(nsample,nvar);
    %method of Stein
    for j=1: nvar
        idx=randperm(nsample);
        P=(idx'-ran(:,j))/nsample;       % probability of the cdf
        s(:,j) = xmean(j) + ltqnorm(P).* xsd(j); % this can be replaced by any inverse distribution function
    end
end


if distrib == 'unif'  % you only need to specify xmin & xmax
    if xmin==0
        xmin=1e-300;
    end
    nvar=length(xmin);
    ran=rand(nsample,nvar);
    s=zeros(nsample,nvar);
    for j=1: nvar
        idx=randperm(nsample);
        P =(idx'-ran(:,j))/nsample;
        xmax(j);
        xmin(j);
        xmax(j)/xmin(j);
        if (xmax(j)<1 & xmin(j)<1) || (xmax(j)>1 & xmin(j)>1)
            'SAME RANGE';
            if (xmax(j)/xmin(j))<threshold %% It uses the log scale if the order of magnitude of [xmax-xmin] is bigger than threshold
                '<1e3: LINEAR SCALE';
                s(:,j) = xmin(j) + P.* (xmax(j)-xmin(j));
            else
                '>=1e3: LOG SCALE';
                s(:,j) = log(xmin(j)) + P.*abs(abs(log(xmax(j)))-abs(log(xmin(j))));
                s(:,j) = exp(s(:,j));
            end
        else
            'e- to e+';
            if (xmax(j)/xmin(j))<threshold %% It uses the log scale if the order of magnitude of [xmax-xmin] is bigger than threshold
                '<1e3: LINEAR SCALE';
                s(:,j) = xmin(j) + P.* (xmax(j)-xmin(j));
            else
                '>=1e3: LOG SCALE';
                s(:,j) = log(xmin(j)) + P.*abs(log(xmax(j))-log(xmin(j)));
                s(:,j) = exp(s(:,j));
            end
        end
    end
end
%hist(s)   % plots the histogram of the pdf

I got this error message

??? Input argument "nsample" is undefined.

Error in ==> LHS_Call at 8

if nsample==1

Please, I need help concerning them.

Answer 1

[Walter Roberson]

You need to run Model_LHS which will call the other routines

Parameter_settings_LHS;

That is a script, and due to changes in the way MATLAB handles functions, when variables are set inside of scripts that are called inside of functions, MATLAB will not always notice that the variable has been set, or might give unusual error messages that have to be read very carefully to see what is happening.

Note: in LHS_Call there are couple of lines that look similar to

if distrib == 'unif'

Notice the single quotes there. single-quoted lists of characters act to create a character vector, which is effectively a uint16 (16 bit integer) numeric array marked as being character oriented. The rules for == for character vectors are the same as for double precision and other numeric types, namely that in order for == to be used, the sizes to compare must be "compatible". For example,

    'abc' == 'b'

asks to compare each element of ['a', 'b', 'c'] to the scalar 'b' returning a 1 x 3 logical vector. And asking

    'abc' == 'abcd'

would be an error because the sizes of the underlying numeric vectors are different in an incompatible way.

== should seldom be used with character vectors when you are comparing more than one character. You should use strcmp() for the comparison instead, such as

  if strcmp(distrib, 'unif')

That version will not give an error if the lengths do not happen to match.

[Chinwendu]

Thank you very @Walter Roberson. It works. I am trying to get PRCC_PLOT and p_value bar plots. I do not know what code I will write in the MODEL_LHS to get them. I am new in using matlab. Thank you

[Walter Roberson]

If you want to change PRCC_PLOT to be a bar plot, then near the bottom of the file where you see the phrase

plot((eval(c4)),(eval(c5)),'.')

in the middle of a line, change that to

bar((eval(c4)),(eval(c5)))

[Torsten]

@Chinwendu comment moved here:

This is the PRCC_PLOT . it shows these error message

>> PRCC_PLOT

Not enough input arguments.

Error in PRCC_PLOT (line 38)

[a k]=size(X); % Define the size of LHS matrix.

I notice that the MODEL_LHS is the the runing script while others you call them. So, the code will be in the MODEL_LHS file. It is how to wirte it that is my propble now.

Thank you very much for your help.

function PRCC_PLOT(X, Y, s, PRCC_var, y_var)
% % % time_points=[200 400]; %why is this hard coded
% % % s=1:length(time_points);
% % % Y=Y(s,:);
[a k]=size(X); % Define the size of LHS matrix
Xranked=rankingN(X);
Yranked=ranking1(Y);
for i=1:k  % Loop for the whole submatrices, Zi
    c1=['LHStemp=Xranked;LHStemp(:,',num2str(i),')=[];Z',num2str(i),'=[ones(a,1) LHStemp];LHStemp=[];'];
    eval(c1);
end
for i=1:k
    c2=['[b',num2str(i),',bint',num2str(i),',r',num2str(i),']= regress(Yranked,Z',num2str(i),');'];
    c3=['[b',num2str(i),',bint',num2str(i),',rx',num2str(i),']= regress(Xranked(:,',num2str(i),'),Z',num2str(i),');'];
    eval(c2);
    eval(c3);
end
for i=1:k
    c4=['r',num2str(i)];
    c5=['rx',num2str(i)];
    [r p]=corr(eval(c4),eval(c5));
    a=['[PRCC , p-value] = ' '[' num2str(r) ' , '  num2str(p) '].'];% ' Time point=' num2str(s-1)];
    figure,bar((eval(c4)),(eval(c5)),'.'),title(a),...
        legend(PRCC_var{i}),xlabel(PRCC_var{i}),ylabel(y_var);%eval(c6);
end

[Chinwendu]

I got another code that plot PRCC_PLOT and P_VALUE but the code shows error messgae too and the plot is coded in MODEL_LHS file. Attached is the file code.

It is how to get the same result that is my propblem. Thank you very much.

[Walter Roberson]

Change

% PRCC_PLOT(LHSmatrix, model.state.(model.analyzeThisOutput).lhs, length(time_points), PRCC_var, model.analyzeThisOutput)

to

PRCC_PLOT(LHSmatrix, model.state.(model.analyzeThisOutput).lhs, length(time_points), PRCC_var, model.analyzeThisOutput)

which is just removing the % that is commenting the line out.

[Chinwendu]

Thank you very much for your reply. I have done that which is in Line 72. But I still have these error messages below. Which I do not understand what is 'indexChosen'.

Thank you for your time.

>> Model_LHS

Unrecognized function or variable 'indexChosen'.

Error in getHeaderContentIndex (line 15)

if(indexChosen == numel(headerIndex))

Error in loadPRCCconfig (line 16)

paramFile = file(getHeaderContentIndex(file,headerStartChar,headerEndChar,'param'));

Error in Model_LHS (line 13)

model = loadPRCCconfig(model,'lhs-prcc-modified/PRCCconfig.txt'); %contains config for PRCC min,baseline,max,initial

[Chinwendu]

This is the code

%% "runs" in the script below)
clear all; close all;
%% [EDITABLE] Sample size N
runs=100;

%% [EDITABLE] Load Model Parameters and Variables
model.analyzeThisOutput = 'V'; % Chosen output variable name to do PRCC analysis

model.dir = pwd; %directory of ODE model and config
model = loadModel(model);
model = loadPRCCconfig(model,'lhs-prcc-modified/PRCCconfig.txt'); %contains config for PRCC min,baseline,max,initial

alpha = 0.001; %0.05; %threshold for significant PRCCs (uncorrelated < alpha)
% Parameter Labels 
PRCC_var = model.paramLabel; %{'s', '\mu_T', 'r', 'k_1','k_2', '\mu_b','N_V', '\mu_V','dummy'};

% Variable Labels
y_var_label = model.yVarLabel; %{'T','T*','T**','V'};

%% [EDITABLE] TIME SPAN OF THE SIMULATION
t_end = 4000; % length of the simulations
tspan = (0:1:t_end);   % time points where the output is calculated
time_points = [2000 4000]; % time points of interest for the US analysis

%% [EDITABLE] INITIAL CONDITION FOR THE ODE MODEL
%  Change the values in the PRCCconfig.txt
% OR set manually: y0(find(strcmp(model.allStateName,'V'))) = 1e-3
y0 = [];
for yi=1:numel(model.allStateName)
    y0 = [y0, model.state.(model.allStateName{yi}).initial]; 
end

for lhsIdx=1:numel(model.paramName)
    model.param.(model.paramName{lhsIdx}).LHS = LHS_Call(model.param.(model.paramName{lhsIdx}).min,...
                                                         model.param.(model.paramName{lhsIdx}).baseline,...
                                                         model.param.(model.paramName{lhsIdx}).max,0, runs, 'unif');
end

%% LHS MATRIX and PARAMETER LABELS
LHSmatrix = [];
for lhsIdx=1:numel(model.paramName)
    LHSmatrix = [LHSmatrix model.param.(model.paramName{lhsIdx}).LHS];
end
for x=1:runs %Run solution x times choosing different values
    f=@ODE_LHS;
    x;
    LHSmatrix(x,:);
    [t,y]=ode15s(@(t,y)f(model,t,y,LHSmatrix,x),tspan,y0,[]); 
    A=[t y]; % [time y]
    %% Save the outputs at ALL time points [tspan]
    %T_lhs(:,x)=Anew(:,1);
    %CD4_lhs(:,x)=Anew(:,2);
    %T1_lhs(:,x)=Anew(:,3);
    %T2_lhs(:,x)=Anew(:,4);
    %V_lhs(:,x)=Anew(:,5);
    
    %% Save only the outputs at the time points of interest [time_points]:
    %% MORE EFFICIENT
    for stIdx=1:numel(model.allStateName)
        model.state.(model.allStateName{stIdx}).lhs(:,x) = A(time_points+1,stIdx+1);
    end
end
%% Save the workspace
% save Model_LHS.mat;
% CALCULATE PRCC
[prcc sign sign_label hFig] = PRCC(LHSmatrix,model,time_points,PRCC_var,alpha,'horizontal'); %'vertical');
% https://www.mathworks.com/matlabcentral/answers/376781-too-many-input-arguments-error

% PRCC_PLOT(X,Y,s,PRCC_var,y_var)
PRCC_PLOT(LHSmatrix, model.state.(model.analyzeThisOutput).lhs, length(time_points), PRCC_var, model.analyzeThisOutput)
prompt = 'Do you want to save all variables and figures? If yes, Enter [Y]. If not enter random: ';
savePrompt = input(prompt, 's');
if(savePrompt == 'Y' || savePrompt == 'y')
    nowVector = clock();
    for tIdx=1:5
        nowCell{tIdx} = num2str(nowVector(tIdx));
    end
    timeString = strjoin(nowCell,'-');
    saveDir = ['results/' 'variable_' model.analyzeThisOutput '/' timeString]; mkdir(saveDir);
    save([saveDir '/PRCCofOutputVar_' model.analyzeThisOutput '.mat'],'-regexp','^(?!(hFig)$).');
    fprintf(['MAT file saved in ' saveDir '\n']);
    saveDir = ['results/' 'variable_' model.analyzeThisOutput '/' timeString '/tables']; mkdir(saveDir);
    for tpIdx=1:numel(time_points)
        fid=fopen([saveDir '/unsortedPRCC' num2str(tpIdx) '_' num2str(time_points(tpIdx)) '.csv'],'w');
        for parIdx=1:numel(model.paramName)
           fprintf(fid,model.paramName{parIdx});
           if(parIdx == numel(model.paramName))
               break;
           else
               fprintf(fid, ',');
           end
        end
        fprintf(fid,'\n');
        for parIdx=1:numel(model.paramName) 
           fprintf(fid,'%.4f',prcc(tpIdx,parIdx)); 
           if(parIdx == numel(model.paramName))
               break;
           else
               fprintf(fid, ',');
           end
        end
        fclose(fid);
    end
    fprintf(['Tables saved in ' saveDir '\n']);
    saveDir = ['results/' 'variable_' model.analyzeThisOutput '/' timeString '/figures'];
    mkdir(saveDir);
    for hfIdx=1:numel(hFig)
        for tpIdx=1:numel(hFig{hfIdx}.figure)
            saveas(hFig{hfIdx}.figure{tpIdx},[saveDir '/' hFig{hfIdx}.name num2str(tpIdx) '_' num2str(time_points(tpIdx)) '.png']);
            saveas(hFig{hfIdx}.figure{tpIdx},[saveDir '/' hFig{hfIdx}.name num2str(tpIdx) '_' num2str(time_points(tpIdx)) '.fig']);
        end
    end
    fprintf(['Figures saved in ' saveDir '\n']);
else
    disp('Figure not saved. If you still want to save, please save them manually.');
end

[Torsten]

Do you really have a file with name "PRCCconfig.txt" in [current directory]/lhs-prcc-modified/ ?

Did you first test the code using the files provided under

http://malthus.micro.med.umich.edu/lab/usadata

?

[Chinwendu Madubueze]

Thank you for your prompt reply. It is from the link. I could not plot PRCC with the code. I have to get another matlab which is from the link below. Thank you very much.

https://github.com/BagaskaraPutra/lhs-prcc-modified

[Walter Roberson]

importdata() is returning something that 'param' cannot be located in.

I never use importdata(), because importdata() can return any of three different datatypes depending on the exact content of the file.

[Chinwendu]

Okay. Thank you very much for your time and reply. i am grateful. I have another code that PRCC_PLOT and it ran but it is not in as bar chart too and it does not produce p_values bar chart. I need help on how to make the plot of the PRCC to be a bar chart and also plot the p_values of the parameters as the plot results in this link https://github.com/BagaskaraPutra/lhs-prcc-modified.

The following is the PRCC and PRCC_PLOT that ran without a bar plot.

The PRCC calling codel

close all
clear variables

file_to_load = 0;

alpha = 0.05;

switch file_to_load
    case 0
        file_title='prcc_Model_LHS'; 
        load([file_title(6:end),'.mat']);
        
end

%V
%[prcc{1},sign{1},sign_label{1}]=PRCC(LHSmatrix,peak_val1,1,PRCC_var,alpha);
%Th1
[prcc{1},sign{1},sign_label{1}]=PRCC(LHSmatrix,peak_val2,1,PRCC_var,alpha);
%Th2
[prcc{2},sign{2},sign_label{2}]=PRCC(LHSmatrix,peak_val3,1,PRCC_var,alpha);
%F
%[prcc{4},sign{4},sign_label{4}]=PRCC(LHSmatrix,peak_val4,1,PRCC_var,alpha);
%I4
%[prcc{5},sign{5},sign_label{5}]=PRCC(LHSmatrix,peak_val5,1,PRCC_var,alpha);
%CTL
%[prcc{6},sign{6},sign_label{6}]=PRCC(LHSmatrix,peak_val6,1,PRCC_var,alpha);
%A
[prcc{3},sign{3},sign_label{3}]=PRCC(LHSmatrix,peak_val7,1,PRCC_var,alpha);
%B
[prcc{4},sign{4},sign_label{4}]=PRCC(LHSmatrix,peak_val8,1,PRCC_var,alpha);


save(file_title,'prcc','sign','sign_label','PRCC_var');

and this is runing PRCC_PLOT code

pos = [2,4.5,6,4.5];
pcols = {1/255*[0,196,255],1/255*[225,12,62],1/255*[55,176,116],1/255*[244,201,107],1/255*[149,52,235],1/255*[176,176,176]};
load('Model_LHS.mat');

            %prclab={'peak_val1','peak_val2','peak_val3','peak_val4','peak_val5','peak_val6','peak_val7','peak_val8'};
            %prctitles={'Viral Load Peak Value','Th_1 Peak Value','Th_2 Peak Value','Interferon Peak Value', 'Interlukin 4 Peak Value', 'CTL Peak Value','Antibody Peak Value','B Cell Peak Value'};
            prclab={'peak_val2','peak_val3','peak_val7','peak_val8 '};
            prctitles={'Th1 Peak Value','Th2 Peak Value','Neutralizing Antibody Peak Value','IgGs+IgM Peak Value '};
            PRCC_var={'\psi', '\eta_v', '\gamma_v', '\eta_{th1}', ...
    '\eta_{th1_f}','\beta_{th1_{i4}}', 'k_{th1}','\gamma_{th1}','\eta_{th2}', '\eta_{th2_f}','\beta_{th2_{f}}',...
    '\gamma_{th2}', ' \eta_{f}', '\gamma_{f}', '\eta_{i4}', '\gamma_{i4}', '\eta_{ctl}', '\eta_{ctl_f}',...
    'k_{ctl}',' \gamma_{ctl}', '\eta_{a}', '\gamma_{a}', '\eta_{b_{th1}}', '\eta_{b_{th2}}', ' \gamma_b '};
                
    
             for k=1:length(prcc)
            % for k=1:3
                fig(k) = figure;
                b=barh(prcc{k}(end,:),'FaceColor','flat');
                b.CData = kron(ones(length(prcc{k}(end,:)),1),pcols{3});
                %b.CData(sign_label{k}.index{end},:)=kron(ones(length(sign_label{k}.index{end}),1),pcols{1});
                %b.CData(sign_label{k}.index{end}(abs(str2num(sign_label{k}.value{end}))>0.5),:)=kron(ones(sum(abs(str2num(sign_label{k}.value{end}))>0.5),1),pcols{1});
                hold on
                %barh(0,'FaceColor',pcols{2});
               % barh(0,'FaceColor',pcols{end});
                hold off
                ax=gca;
                ax.XAxis.LineWidth = 3;
                ax.YAxis.LineWidth = 3;
                %ax.FontSize = 20;
                yticks([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])
                yticklabels(PRCC_var);
                set(gca,'fontweight','bold','fontsize',15) 
                fontSize1 = 30;
                fontSize2 = 30;
                xlim([-1.01,1.01]);
                figprop(k).xlab=xlabel('PRCC','FontSize', fontSize1,'fontname','times');
                %figprop(k).leg = legend({strcat('$p<0.05$,', 10, '$\rm{PRCC}<0.5$'),strcat('$p<0.05$,',10,'$\rm{PRCC}>0.5$'),'$p>0.05$'},'location','best');
                figprop(k).ax = gca;
                figprop(k).title=title(prctitles{k},'FontSize', fontSize2,'fontname','times');
                figprop(k).fnt_size=20;
                %f = figure;  
                figprop(k).Position = [10 10 550 4000];
                %saveas(fig(k),'Desktop\fig(k).pdf')

            end

this is one of the figure the code above produce.

but I want to get something like the following plots

How do I get P_values plot of the parameters from the code. Why is the code not producing a bar plot like this above?

Thank you very much for answering me. I am grateful.