%Code will produce a plot showing the number of infants required to achieve a 
% power of 95% for different noise levels for a given intervention effect 
% comparing 2 groups, with and without an intervention. 
% Here the term noise level is used to refer to the standard deviation of
% the residuals in the correlation between nociceptive sensitivity and the
% response to the painful procedure.
%The power is computed when the groups are compared with and without 
% taking into account individual nociceptive sensitivity. 
% Simulations are run with increasing sample size until a power of 95% is achieved.
% As Figure 2E.

% Please note that this code takes some time to run if the noise level is high. 
% To replicate the figures in the paper exactly please use [0.5:0.5:5].
% The level of noise can be adjusted on line 33
% Progress on the simulations is displayed in the command window.
% All simulations with nociceptive senstivity are performed first, followed
% by those without accounting for nociceptive senstivity. 

% The code uses an intervention effect of 40% reduction. To change this
% please alter the intervention_effect variable on line 32.

% Code produced by Caroline Hartley August 2020
% Please see the related paper and cite as:
%


%%


intervention_effect=40; %percentage intervention effect
noise_level=[0.5:0.5:5]; %note the noise level is a multiple of 0.37 - the standard deviation of residuals in Study 1

%initialise variables to store the number of infants needed per noise level
% to reach a power of 95%
number_for_95power_with_nociceptivesensitivity=zeros(size(noise_level));
number_for_95power_no_nociceptivesensitivity=zeros(size(noise_level));

drug_reduction=intervention_effect;
iter_n=1000;

%run loop taking into account nociceptive senstivity
% for each noise level the number of subjects per group is increased by 1
% (starting with one infant per group) until 95 power is reached

for m=1:length(noise_level)
    disp(strcat('running simulations with nociceptive sensitivity and noise level of: ',num2str(noise_level(m))))
    power=0; %initialise power variable
    no_subs=1;
    
    while power<95
        no_subs=no_subs+1;
        n_p=no_subs; n_a=n_p;

        noise_level_iter=noise_level(m);

        p_store_with_nociceptivesensitivity=zeros(iter_n,1);

        for iter=1:iter_n
            rand('seed',iter) %seed variables for replicability
            randn('seed',iter)

            %simulate individual nociceptive sensitivity (pinprick responses)
            pp_p=1.4*rand(n_p,1)+0.15; %placebo data simulated nociceptive sensitivity
            pp_a=1.4*rand(n_a,1)+0.15; %analgesic data simulated nociceptive sensitivity

            %simulate lance responses from these nociceptive sensitivty values
            %Related to the relationship found in Study 1 - see Methods
            l_p_estimate=2.62*pp_p-0.75; %placebo lance data - estimated data from linear regression
            l_a_estimate=2.62*pp_a-0.75; %analgesic lance data

            %reduce l_a data by drug effect
            l_a=(1-(drug_reduction/100))*l_a_estimate;

            %add noise
            noise=noise_level_iter*0.37; %noise relates to standard deviation of residuals. Standard deviation of residuals in study 1 was 0.37
            l_p=l_p_estimate+noise*randn(n_p,1);
            l_a=l_a+noise*randn(n_a,1);

            %compare data accounting for individual nociceptive sensitivity
            lance=[l_p;l_a];
            pp=[pp_p;pp_a];
            modality=[ones(length(l_p),1);2*ones(length(l_a),1)];
            tbl=table(lance,pp,modality,'VariableNames',{'Lance','Pinprick','modality'});
            tbl.modality = categorical(tbl.modality);
            lm = fitlm(tbl,'Lance~Pinprick+modality');
            
            p_store_with_nociceptivesensitivity(iter)=table2array(lm.Coefficients(3,4));


        end

        %find the power
        power=length(find(p_store_with_nociceptivesensitivity<0.05))/iter_n*100; 

    end

	number_for_95power_with_nociceptivesensitivity(m)=no_subs;  
    
end


%run loop without taking into account nociceptive senstivity
% for each noise level the number of subjects per group is increased by 1
% (starting with one infant per group) until 95 power is reached

for m=1:length(noise_level)
    disp(strcat('running simulations without nociceptive sensitivity and noise level of: ',num2str(noise_level(m))))
    power=0; %initialise power variable
    no_subs=1;
    
    while power<95
        no_subs=no_subs+1;
        n_p=no_subs; n_a=n_p;

        noise_level_iter=noise_level(m);

        p_store_no_nociceptivesensitivity=zeros(iter_n,1);

        for iter=1:iter_n
            rand('seed',iter) %seed variables for replicability
            randn('seed',iter)

            %simulate individual nociceptive sensitivity (pinprick responses)
            pp_p=1.4*rand(n_p,1)+0.15; %placebo data simulated nociceptive sensitivity
            pp_a=1.4*rand(n_a,1)+0.15; %analgesic data simulated nociceptive sensitivity

            %simulate lance responses from these nociceptive sensitivty values
            %Related to the relationship found in Study 1 - see Methods
            l_p_estimate=2.62*pp_p-0.75; %placebo lance data - estimated data from linear regression
            l_a_estimate=2.62*pp_a-0.75; %analgesic lance data

            %reduce l_a data by drug effect
            l_a=(1-(drug_reduction/100))*l_a_estimate;

            %add noise
            noise=noise_level_iter*0.37; %noise relates to standard deviation of residuals. Standard deviation of residuals in study 1 was 0.37
            l_p=l_p_estimate+noise*randn(n_p,1);
            l_a=l_a+noise*randn(n_a,1);

           %compare data without accounting for individual nociceptive sensitivity
            [~,p]=ttest2(l_p,l_a);
            
            p_store_no_nociceptivesensitivity(iter)=p;


        end

        %find the power
        power=length(find(p_store_no_nociceptivesensitivity<0.05))/iter_n*100; 

    end

	number_for_95power_no_nociceptivesensitivity(m)=no_subs;  
    
end


%%
approx_sd_resid=noise_level*0.37; %sd related to sd of data from Study 1

figure; plot(approx_sd_resid,number_for_95power_no_nociceptivesensitivity,'b')
hold on
plot(approx_sd_resid,number_for_95power_with_nociceptivesensitivity,'r')
ylabel('Number of infants','fontsize',15)
xlabel('Standard deviation of residuals','fontsize',15)
set(gca,'fontsize',15)
xlim([min(approx_sd_resid),max(approx_sd_resid)])
legend('Without nociceptive sensitivity','With nociceptive sensitivity')

figure; plot(approx_sd_resid,(number_for_95power_no_nociceptivesensitivity-number_for_95power_with_nociceptivesensitivity)./number_for_95power_no_nociceptivesensitivity*100,'k')
hold on
ylabel('Percentage reduction','fontsize',15)
xlabel('Standard deviation of residuals','fontsize',15)
set(gca,'fontsize',15)
xlim([min(approx_sd_resid),max(approx_sd_resid)])
set(gca,'XTick',approx_sd_resid)