clear all
close all
% Hemopexin in DMD Cases I. Refer to data set dmd.dat|mat|xls from
% Exercise 2.19. The measurements of hemopexin are assumed normal.
%
% (a) Form a 95% confidence interval for the mean response of hemopexin
% h in a population of all female DMD carriers (carrier=1).
% Although the level of pyruvate kinase seems to be the strongest single
% predictor of DMD, it is an expensive measure. Instead, we will explore
% the level of hemopexin, a protein that protects the body from oxidative
% damage. The level of hemopexin in a general population of women of
% the same age range as that in the study is believed to be 85 on average.
%
% (b) Test the hypothesis that the mean level of hemopexin in the population
% of women DMD carriers significantly exceeds 85. Use alpha = 5%.
% Report the p-value as well.
%
% (c)What is the power of the test in (b) against the alternative H1 : mu1 =89.
%
% (d) The data for this exercise come from a study conducted in Canada.
% If you wanted to replicate the test in the USA, what sample size would
% guarantee a power of 99% if H0 were to be rejected whenever the difference
% from the true mean was 4, (|mu0 - mu1| = 4)? A small pilot study
% conducted to assess the variability of hemopexin level estimated the standard
% deviation as s = 12.
load 'C:\STAT\Descriptive\Descriptivedat\dmd.mat'
alpha = 0.05;
hemo = dmd( dmd(:,6) == 1, 3);
n=length(hemo)  %n=67
xbar = mean(hemo) %xbar = 92.9303
s2  = var(hemo) %97.1721
s = sqrt(s2)  %9.8576
mu0 = 85;
t = (xbar-mu0)/sqrt(s2/n)  %6.5850


%(a) 
[xbar - tinv(1-alpha/2,n-1) * s/sqrt(n),... 
    xbar + tinv(1-alpha/2,n-1) * s/sqrt(n)] 
%[90.5258   95.3347]

%(b)
tcrit = tinv(1-alpha,n-1) %1.6683
 
pval = 1-tcdf(t, n-1) 
%pval =  4.4003e-009

%(c)
mu1 =89;
pow1 = @(n) nctcdf( -tinv(1-alpha, n-1), n-1,-abs(mu0-mu1)*sqrt(n)/s);
pow1(n)
% 0.9497

beta = 0.01;
pow2 = @(n) nctcdf( -tinv(1-alpha, n-1), n-1,-abs(mu0-mu1)*sqrt(n)/s) - (1-beta)
%(d)

approxssize =  s^2 *(norminv(1-alpha) + norminv(1-beta))^2/(mu0 - mu1)^2  %95.7779
exactssize = fzero(pow2, approxssize)  %97.1509





