%Code will produce a plot of the power at different sample sizes for
%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. As Figure 2B.

% 1000 samples (line 37) are compared for each sample size, and the power calculated
% as the percentage of significant (p<0.05) results at each sample size.

% The level of intervention effect can be adjusted on line 16
% The number of subjects which are included can be adjusted on line 17

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

drug_reduction=40; %percentage drug reduction or other intervention effect
no_subs=[5:10:80]; %number of subjects per group

%% Simulate


%test if lance results are different between placebo and anagesic groups by
%themselves and whether different if take into account sensitivty

%assign variables to store the power for each sample size
power_no_nociceptivesensitivity=zeros(1,length(no_subs));
power_with_nociceptivesensitivity=zeros(1,length(no_subs));


for j=1:length(no_subs) %loop through all sample sizes

%simulate with n_p subjects with 'placebo' and 'n_a' subjects with
%'analgesic'/intervention

n_p=no_subs(j); n_a=n_p;

noise_level_iter=1; %this relates to the noise level. Setting it to 1 means that the standard deviation of residuals will be the same as Study 1

iter_n=1000; %1000 simulated data sets for each sample size

%set-up variables to store info for all iterations - will store p-values

p_store_no_nociceptivesensitivity=zeros(iter_n,1);
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 without accounting for individual nociceptive sensitivity
[~,p]=ttest2(l_p,l_a);
p_store_no_nociceptivesensitivity(iter)=p;

%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

power_no_nociceptivesensitivity(j)=length(find(p_store_no_nociceptivesensitivity<0.05))/iter_n*100;
power_with_nociceptivesensitivity(j)=length(find(p_store_with_nociceptivesensitivity<0.05))/iter_n*100;

end

%%
figure; plot(no_subs,power_no_nociceptivesensitivity,'b')
hold on
plot(no_subs,power_with_nociceptivesensitivity,'r')
legend('Without nociceptive sensitivity','With nociceptive sensitivity')
xlabel('Number of infants per group','fontsize',15)
ylabel('Power','fontsize',15)
set(gca,'fontsize',15)