%---------------MFILE ARTHRITIS: Comploglog Link---------------------
% Arthritis Treatment Data. The data were obtained from Koch and Edwards (1988)
% for a double-blind clinical trial investigating a new treatment for rheumatoid
% arthritis. In this data set, there were 84 subjects with different ages who
% received an active or placebo treatment for their arthritis pain, and the
% subsequent extent of improvement was recorded as marked, some, or none.
% The dependent variable "improve" was an ordinal categorical observation
% with three categories (0 for none, 1 for some, and 2 for marked).
% The three explanatory variables were "treatment" (1 for active or 0
% for placebo), "gender" (1 for male, 2 for female), and "age"
% (recorded as a continuous variable). The data in arthritis2.dat is
% organized as a matrix, with 84 rows corresponding to subjects and 5 columns
% containing id number, treatment, gender, age, and improvement status.
%
% 57 1 1 27 1
% 9 0 1 37 0
% 46 1 1 29 0
% ...
% 15 0 2 66 1
% 1 0 2 74 2
% 71 0 2 68 1
%
% Dichotomize the variable "improve" as "improve01=improve>0;"
% Fit the binary regression with "improve01" as a response and "treatment",
% "gender", and "age" as covariates. Use the three links: logit,
% probit and comploglog and compare the models by comparing deviances
% and model residuals.
clear all
close all
disp('Binary Regression: Arthritis2')
lw = 2.5;
set(0, 'DefaultAxesFontSize', 16);
fs = 15;
msize = 10;
%-------------------------------------------------------------------------
load 'C:\STAT\Logistic\Logisticdat\arthritis2.dat'
caseid = arthritis2(:,1);
treatment = arthritis2(:,2);
gender = arthritis2(:,3);
age = arthritis2(:,4);
improve = arthritis2(:,5);
improve01 = arthritis2(:,5)>0 ;
X = [treatment gender age];
[betas, deviance, stats]=glmfit(X,improve01,'binomial','link','comploglog')
figure(1)
score = betas(1) + betas(2)*treatment + betas(3)*gender + betas(4)* age;
plot(score, improve01,'o',...
'MarkerSize',msize, 'MarkerEdgeColor','k', 'MarkerFaceColor','g')
xx = -2.7:0.01:1.4;
imp = 1 - exp(- exp(xx) );
hold on
plot(xx, imp,'r-','LineWidth',lw)
xlabel('Score')
ylabel('Probability of Improving')
%%
[betas2, deviance2, stats2]=glmfit(X,improve01,'binomial','link','logit')
[betas3, deviance3, stats3]=glmfit(X,improve01,'binomial','link','probit')
score2 = betas2(1) + betas2(2)*treatment + betas2(3)*gender + betas2(4)* age;
score3 = betas3(1) + betas3(2)*treatment + betas3(3)*gender + betas3(4)* age;
figure(2)
imp = 1 - exp(- exp(score) );
imp2 = exp(score2)./(1 + exp(score2));
imp3 = normcdf(score3);
plot(score, imp,'r*')
hold on
plot(score2, imp2,'ko')
plot(score3, imp3,'bd')
xlabel('scores')
ylabel('fits')
legend('cloglog','logit','probit',2)
deviance %92.0751
deviance2 %92.0628
deviance3 %91.9286