%Code will produce a plot showing the number of infants required to achieve a 
% power of 95% for different levels of intervention effect comparing 2 groups, 
% with and without an intervention. 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 2C and D.

% Please note that this code takes some time to run if the intervention
% effect is low. To replicate the figures in the paper exactly please use [5:5:95]
% The level of intervention effect can be adjusted on line 24.
% 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. 

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


%%

intervention_effect=[5:5:95]; %percentage intervention effect

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


%run loop taking into account nociceptive senstivity
% for each drug 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(intervention_effect)
    disp(strcat('running simulations with nociceptive sensitivity and intervention effect of: ',num2str(intervention_effect(m))))
    power=0; %initialise power variable
    
    no_subs=1; %start with 1 infant per group
    
    while power<95
        
        no_subs=no_subs+1; %increase number of subjects by 1 each time until power is 95
        n_p=no_subs; n_a=n_p;

        noise_level_iter=1;
        drug_reduction=intervention_effect(m);
        iter_n=1000;
        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 taking into account nociceptive senstivity
% for each drug 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(intervention_effect)
    disp(strcat('running simulations without nociceptive sensitivity and intervention effect of: ',num2str(intervention_effect(m))))
    power=0; %initialise power variable
    
    no_subs=1; %start with 1 infant per group
    
    while power<95
        
        no_subs=no_subs+1; %increase number of subjects by 1 each time until power is 95
        n_p=no_subs; n_a=n_p;

        noise_level_iter=1;
        drug_reduction=intervention_effect(m);
        iter_n=1000;
        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

%%
figure; semilogy(intervention_effect,number_for_95power_no_nociceptivesensitivity,'b')
hold on
semilogy(intervention_effect,number_for_95power_with_nociceptivesensitivity,'r')
ylabel('Number of infants','fontsize',15)
xlabel('Intervention effect','fontsize',15)
set(gca,'fontsize',15)
xlim([min(intervention_effect),max(intervention_effect)])
legend('Without nociceptive sensitivity','With nociceptive sensitivity')

figure; plot(intervention_effect,(number_for_95power_no_nociceptivesensitivity-number_for_95power_with_nociceptivesensitivity)./number_for_95power_no_nociceptivesensitivity*100,'k')
ylabel('Percentage reduction','fontsize',15)
xlabel('Intervention effect','fontsize',15)
set(gca,'fontsize',15)
xlim([min(intervention_effect),max(intervention_effect)])