##### https://gitlab.com/paediatric_neuroimaging/simulating_power_nociceptive_sensitivity

Tip revision:

**1a465a40228c1c5ab5d33d4cbdf8cd99b29e9fcf**authored by**Maria M Cobo**on**21 April 2021, 14:31 UTC****Update README.md** Tip revision:

**1a465a4** calculate_number_of_subjects_for_power95_with_diff_interaction.m

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