%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)