 
%Second birth CPH
close all force
clear all 
disp('Time Between 1st and 2nd Birth')
  lw = 2; 
  set(0, 'DefaultAxesFontSize', 16);
  fs = 15;
  msize = 10;
% end of preabmle--------------------------------------
%
% The Medical Birth Registry of Norway was established in 1967 and contains 
% information on all births in Norway since that time. 
% The data set   secondbirth.dat|xlsx  distilled from the registry provides 
% the time between first and second births for a selection of 53,558 women. 
% It is of  interest to explore whether this time was possibly affected if the
% first child died within one year of its birth.
% 
% The data set contains the following variables (as columns)
% 
%  age --   age of mother at first birth (in years)  
%  sex --    sex of first child (1 = boy, 2 = girl) 
%  death --  first child died within one year (0 = no, 1  =yes)
%  time  --   time from first birth to second birth or censoring (in days) 
%  status --    censoring indicator (0 = censored, 1 = birth)  
% 
 % 
% (a) Fit the Cox proportional hazard regression
% with variables  age,  sex, and  death  as covariates.
% What is the 95\% CI for the slope  beta_3, corresponding to 
% variable   death ?
% Is the variable   sex  significant? 
% 
% (b) Plot the hazard function (alternatively, cummulative
%  hazard or survival function) for    death = 0  and
%  death = 1 .
%  
% (c) Using Mantel's procedure  manteltwo.m, test the hypothesis that
% the two survival functions are the same.
% 
% 
 %  Reference: The Medical Birth Registry of Norway is acknowledged for 
%allowing the usage of the data and Dr. Stein Emil Vollset for providing the data.
 
%(a)
a = importdata('../Survivaldat/secondbirth.dat')  
age=a.data(:,1);
sex =a.data(:,2);
death=a.data(:,3);
status=a.data(:,5);
censor=1-status;
%
time = a.data(:,4);
% 
ple(time(death==0),censor(death==0))
hold on
ple(time(death==1),censor(death==1))
hold off
% 
    [b,logl,H,stats] = coxphfit([age sex death],time,...
    'censoring',censor);
%
% stats.covb
%     0.0000   -0.0000    0.0000
%    -0.0000    0.0002    0.0000
%     0.0000    0.0000    0.0066
% 
% stats.beta
%     0.0877
%    -0.0187
%     0.8889
% 
% stats.se
%    0.0036
%     0.0157
%     0.0813
% 
% stats.z
%    24.3223
%    -1.1906;
%    10.9364
% 
% stats.p
%     0.0000
%     0.2338
%     0.0000
% 
 [stats.beta(3) - norminv(0.975)*stats.se(3)  stats.beta(3) + norminv(0.975)*stats.se(3)]
%   0.7296    1.0482

%(b)
figure(1)
stairs(H(:,1),H(:,2))
hold on
stairs(H(:,1),exp(0.8889)*H(:,2),'r-')
legend('Baseline','Death 1st Child','Location','NorthWest' )
%or
figure(2)
stairs(H(1:length(H(:,1))-1,1),diff(H(:,2)))
hold on
stairs(H(1:length(H(:,1))-1,1),diff(exp(0.8889)*H(:,2)),'r-')
legend('Baseline','Death 1st Child','Location','NorthWest' )
%or
figure(3)
stairs(H(:,1),exp(-H(:,2)))
hold on
stairs(H(:,1),exp(-exp(0.8889)*H(:,2)),'r-')
legend('Baseline','Death 1st Child')

% 
%(c)
tA= time(death ==0);
tB= time(death ==1);

cA= censor(death ==0);
cB= censor(death ==1);
      [pval chi2]=manteltwo(tA,  cA, tB, cB)

% pval =
% 
%      0
% 
% 
% chi2 =
% 
%   113.1702



