Poisson vs. Logistic regression

data melanoma;
      input age $ region $ cases total;
      ltotal=log(total);
      datalines;
   35-44 south 75  220407
   45-54 south 68  198119
   55-64 south 63  134084
   65-74 south 45   70708
   75+   south 27   34233
   <35   south 64 1074246
   35-44 north 76  564535
   45-54 north 98  592983
   55-64 north 104 450740
   65-74 north 63  270908
   75+   north 80  161850
   <35   north 61 2880262
   ;
proc genmod data=melanoma order=data;
  class age region;
  model cases = age region / dist=poisson link=log offset=ltotal;
run;

proc sql;
  create table melanoma2 as
  select age, region, 1 as resp, cases as count  from melanoma union
  select age, region, 0 as resp, total-cases as count  from melanoma 
;quit;

proc genmod data=melanoma2;
  weight count;
  class age region;
  model resp=age region / dist=binomial link=logit;
run;

proc mixed

data bond;
  input ingot metal $ pres @@;
  datalines;
1   n   67.0  1   i   71.9  1   c   72.2
2   n   67.5  2   i   68.8  2   c   66.4
3   n   76.0  3   i   82.6  3   c   74.5
4   n   72.7  4   i   78.1  4   c   67.3
5   n   73.1  5   i   74.2  5   c   73.2
6   n   65.8  6   i   70.8  6   c   68.7
7   n   75.6  7   i   84.9  7   c   69.0
;

proc mixed data=bond method=reml;
   class ingot metal;
   model pres=metal;
   random ingot;
   lsmeans metal / diff=control('n') adjust=dunnett;
run;

/* by simulation */
lsmeans metal / diff=control("n") cl
                adjust=simulate(report seed=4943838 cvadjust);

format

/* formats.sas7bcat is saved under h:\temp */

libname ex 'h:\temp';

proc format library=ex;
  value sev 1='mile'
            2='moderate'
            3='severe'; 
run;

/* To use the permanent format */
libname ff 'h:\temp';
options fmtsearch=(ff);

SAS/IML - something worth to remember

# sasdataset -> matrix
use sas_data;
read all into M; /* read all variables into a matrix M;
print M;

# matrix -> sasdataset
proc iml;
M={1 2 3, 4 5 6, 7 8 9};

create sas_data from M;
append from M
close sas_data

quit;

# pi
pi=constant('PI');

# rank of a matrix
rank=round(trace(ginv(a)*a));

proc glmpower

data exemplary;
  do variety=1 to 2;
    do exposure=1 to 3;
      input height @@; output;
    end;
  end;
datalines; /* expected population means */
14 16 21 
10 15 16
run;

proc glmpower data=exemplary;
  class variety exposure;
  model height=variety|exposure;
  contrast "variety" variety 1 -1;
  contrast "exp 1 vs 3" exposure 1 0 -1;
  contrast "inter" variety*exposure 1 1 1, -1 -1 -1;
  power stddev=5
  ntotal=60
  power=.;
  plot x=n min=30 max=90;
run;

proc power - sample size & power

/* one sample t-test */

* power = ? ;
proc power;
onesamplemeans
alpha=0.05
sides=2
nullm=20
mean=22
stddev=4
ntotal=44
power=.;
run;

* sample size = ? ;
proc power;
onesamplemeans
alpha=0.01
sides=u /* U: upper one-sided, L: lower one-sided
nullm=20
mean=22
stddev=4
ntotal=.
power=0.8;
run;

/* paired t-test */
proc power;
pairedmeans test=diff
alpha=.01
sides=2
meandiff=3
stddev=3.5
corr=.2
npaires=20 30 40
power=.;
run;

/* independent t-test */
proc power;
twosamplemeans
meandiff=3 to 4 by .5 /* the same as 3 3.5 4 */
stddev=8 to 9 by .5
groupweights=(1 1)
power=0.8
ntotal=.;
plot y=power min=0.5 max=0.99;
run;

* a different way;
proc power;
twosamplemeans
groupmeans=(13 14) (13 14.5) (13 15) /* same as 13|14 14.5 15 */
stddev=1.2 1.7
groupweights=1|1 2 3 /* same as (1 1) (1 2) (1 3) */
power=0.8
ntotal=.;
run;

* Power vs Effect Size
proc power;
twosamplemeans test=diff
meandiff=0 to 2.5 by 0.5
stddev=.5657 1.0 1.4318
power=.
npergroup=10;
plot x=effect interpol=join;
run;

/* Multiple Regression */
proc power;
multreg
model=random
nfullpredictors=7
ntestpredictors=1
partialcorr=0.35
ntotal=100
power=.;
plot x=n min=50 max=150;
run;

/* One-way ANOVA */
proc power;
onewayanova test=overall
alpha=.05
groupmeans=(5 7 3 11)
stddev=4 5 6
npergroup=10 15
power=.;
run;

/* Normal Approximation to test a proportion */
proc power;
onesamplefreq test=z method=normal /* test=adjz with continuity corrrection */
sides=1 /* one-sided */
alpha=.05
nullproportion=0.3
proportion=.2
ntotal=.
power=.8;
run;

/* Fisher's exact test */
proc power;
twosamplefreq test=fisher
proportiondiff=0.10 to 0.15 by 0.01
refproportion=.2
npergroup=150
power=.;
run;

/* LR Chi-square Test for Two Proportions */
/* test=pchi for Pearson Chi-square Test for Two Proportions */
proc power;
twosamplefreq test=lrchi
proportiondiff=0.10 to 0.15 by 0.01
refproportion=.2
npergroup=150
power=.;
run;

/* Correlation */
proc power;
onecorr dist=fisherz
npvars=6
corr=.35
nullcorr=.2
sides=1
ntotal=100
power=.;
run;

/* comparing 2 survival curves */
proc power;
twosamplesurvival test=logrank
gexphs=0.3567 | 0.5978 .6931
grouplossexphazards=(0.3567 0.3567)
accrualtime=1
followuptime=1
groupweights=(1 2)
power=.
ntotal=225;
run;

/* TOST */
proc power;
twosamplemeans test=equiv_ratio
lower=.8
upper=1.25
meanratio=1 1.2
cv=.1 .2 .3
npergroup=.
power= .8 .9;
run;

Simple ODS

# Trace output
ods trace on / label;

proc mixed;
...
run;

ods trace off;

# use it
ods select none;
ods output LSMeans=lsm;
proc mixed;
...
run;
ods select all;