%tpower2.m
%create t and noncentral t distribution figures to be used in power lecture

clear all;

%control group stats
mean_control = 1018;
sd_control = 192;
n_control = 20;
%treatment group
mean_treatment = 880;
sd_treatment = 141;
n_treatment = n_control;

%conventional alpha level
alpha = 0.05;
%degrees of freedom (t)
df = n_treatment+n_control-2;
%effect size (Cohen's 'd' as defined in G*Power)
d = (mean_treatment-mean_control)/((sd_treatment+sd_control)/2);
%non-centrality parameter
delta = d * sqrt((n_treatment*n_control)/(n_treatment+n_control));

%critical t-value 
tcrit = tinv(alpha, df);
%derive beta and power probabilities (for treatment < control)
ePower = nctcdf(tcrit, df, delta);
beta = 1-ePower;
 
%plot central and noncentral t distributions
t = [-8:0.1:8];
tdist = tpdf(t,df);
nct = nctpdf(t,df,delta);
%central t distribution
figure(1), clf, plot(t,tdist,'b-','LineWidth',2);
axis([-8 8 0 0.5]);
xlabel('t-value'), ylabel('Probability Density');
title(['t Distribution  (df=',num2str(df),')']);

%central t with t-critical value
figure(2), clf, plot(t,tdist,'b-','LineWidth',2);
xlabel('t-value'), ylabel('Probability Density');
title(['t Distribution  (df=',num2str(df),')']);
hold on;
plot(t,tdist,'b-','LineWidth',2);
line([tcrit tcrit],[0 0.5],'LineWidth',2,'Color',[0 0 0]);

%noncentral t
figure(3), clf, plot(t,nct,'g-','LineWidth',2);
xlabel('t-value'), ylabel('Probability Density');
title(['noncentral t(38,-2.6)']);
hold on;
line([tcrit tcrit],[0 0.5],'LineWidth',2,'Color',[0 0 0]);

%central and non-central t distributions
figure(4), clf, plot(t,nct,'g-','LineWidth',2);
xlabel('t-value'), ylabel('Probability Density');
title(['noncentral t(38,-2.6) versus t(38)']);
hold on;
plot(t,tdist,'b-','LineWidth',2);
line([tcrit tcrit],[0 0.5],'LineWidth',2,'Color',[0 0 0]);
plot(t,nct,'g-','LineWidth',2);
legend('nct','t');

%annotated central and non-central t distributions
figure(5), clf, plot(t,nct,'g-','LineWidth',2);
xlabel('t-value'), ylabel('Probability Density');
title(['noncentral t(38,-2.6) versus t(38)']);
hold on;
plot(t,tdist,'b-','LineWidth',2);
line([tcrit tcrit],[0 0.5],'LineWidth',2,'Color',[0 0 0]);
plot(t,nct,'g-','LineWidth',2);
legend('nct','t');
%annotate figure
fstr(1)={['t-critical=',num2str(tcrit),' (alpha=',num2str(alpha),')']};
fstr(2)={['df=',num2str(df),',  nc-parameter=',num2str(delta)]};
fstr(3)={['Effect Size (d) = ',num2str(d)]};
fstr(4)={['Beta = ',num2str(beta)]};
fstr(5)={['Observed Power = ',num2str(ePower)]};
text(1.6,0.25,fstr);
hold off;
