#***********************************************************
#
# Macro:             BCtrans.MAC
# Version:           Release 14
# Final Date:        04/02/04
# Authors:           Steven Orlich and Mike Delozier
#                    Minitab Inc.
#                    Quality Plaza
#                    1829 Pine Hall Road 
#                    State College, PA  16801-3008
#***********************************************************

#***********************************************************
#
# Neither Minitab Inc. nor the authors of this MACRO makes
# any claim or offers any warranty whatsoever with regard to
# the accuracy of this MACRO or its suitability for use, and
# Minitab Inc. and the authors of this MACRO each disclaims
# any liability with respect thereto.
#
#***********************************************************

#***************************************************************************
#           BCtrans
#       Command and Subcommands

#1.  BCtrans C C...C: Enter the response variable column first followed
#                       by the predictor variable columns.
#2.  range k k:       Use to specify the range for lambda in the PRESS
#                       plot. This range for lambda overrides the default
#                       range that covers the 95% confidence interval.
#3.  influence:       Use to display an index plot showing the influence
#                       of individual cases on the likelihood estimate of
#                       the response transformation parameter.
#4.  bcstore C C:     Use to store the log-likelihood and corresponding
#                       lambda values used in creating the log-likelihood
#                       plot.
#5.  infstore C C:    Use to store the approximate likelihood distances
#                       and case numbers used in creating the index plot.
#6.  presstore C C:   Use to store PRESS and corresponding lambda values
#                       used in creating the PRESS plot.

#***************************************************************************

Macro
BCtrans Y X.1-X.n;
range r1 r2;
influence;
bcstore St1 St2;
infstore St3 St4;
presstore St5 St6.

Mcolumns Y X.1-X.n Y1 Lambda Res Trans LL Temp1 Temp2 Junk
Mcolumns Conf Col Col1 Col2 Col3 Col4 Col5 Final Case Mesh
Mcolumns St1 St2 St3 St4 St5 St6 Oper YY XX.1-XX.n Lev Press
Mcolumns HH.1-HH.100 II.1 Bug JJ.1 T1 T2

Mconstants a c e f g h i j k l m t u lo hi lam err one two
Mconstants aa r1 r2 ss tt ww zz bg

Mmatrix m1 m2 m3 m4 m5 m6 m7 m8 m9

noecho
brief
note
note Macro is running ... please wait
mtitle "Box-Cox Power Transformation Analysis"
note
notitle
brief 0

#***************************************************************************
#   Error checks

min Y err
IF err <= 0
  brief 2
  Note **Error** Response values must be positive.
  Note
  exit
ENDIF

IF infstore and not influence
  brief 2
  Note **Error** Must select both influence and infstore for storage.
  Note
  exit
ENDIF

#***************************************************************************
#   Check for missing values

RNmissing Y X.1-X.n Oper
copy Y X.1-X.n YY XX.1-XX.n;
  omit Oper = 1:1000.             #1000 was chosen as an extreme upper limit

IF max(Oper) > 0
  brief 2
  Note **Note** Rows with missing values have been deleted from the analysis.
  Note
  brief 0
ENDIF

#***************************************************************************
#   Session window information

kkname ss Y
kkset tt "   "
kkcat tt ss tt
name tt 'Response:'

DO i = 1:n
    IF i < n
      kkname ww X.i
      kkset zz " , "
      kkcat ww zz ww
      copy ww HH.i
    ELSE
      kkname ww X.i
      copy ww HH.i
    ENDIF
ENDDO

let II.1[1] = 1
code (1) "Predictor(s):   " II.1 II.1

IF n > 5
  kkset bg "..."                #abbreviate names of predictors if n > 5
  copy bg Bug
  conc II.1 HH.1-HH.4 Bug JJ.1
ELSE
  conc II.1 HH.1-HH.n JJ.1
ENDIF

brief
Note Model Information
Note ------------------------------
Note
print tt
print JJ.1;
format(A80).
Note ------------------------------
brief 0

#***************************************************************************
#   Find optimal lambda based on log-likelihood

loge YY Y1                  #LOGE(Response variable)
set Lambda                  #Lambda Values
-2:2/0.01
end

DO i = 1:401
  let c = Lambda[i]
    IF c = 0
      regress Y1 n XX.1-XX.n;
      residuals Res;
          xpxinv m1.
    ELSE
      let Trans = ((YY**c)-1)/c
      regress Trans n XX.1-XX.n;
        residuals Res;
        xpxinv m1.
    ENDIF
  let LL[i] = ((c-1)*sum(Y1))-(0.5*(count(YY))*(loge((ssq(Res))/(count(YY)))))
ENDDO

sort LL Lambda Temp1 Temp2
max Temp1 e
let f = e-1.92
code (f:e)1 (-9999:f) '*' Temp1 Conf
let Conf = Conf*Temp2
let lam = Temp2[401]
min Conf lo
max Conf hi

#***************************************************************************
#   Compute PRESS for refined lambda range

IF range
  copy r1 one
  copy r2 two
ELSE
  let one = round(lo,1)
  let two = round(hi,1)
    IF lo <= one
      let one = one-0.10
    ENDIF

    IF hi > two
      let two = two+0.10
    ENDIF
ENDIF

set Mesh
one:two/.01
end

count Mesh aa
DO i = 1:aa
  let c = Mesh[i]
    IF c = 0
      regress Y1 n XX.1-XX.n;
        residuals Res;
        hi Lev.
      let Press[i] = ssq(YY-(expo(Y1-(Res/(1-Lev)))))
    ELSE
      let Trans = YY**c
      regress Trans n XX.1-XX.n;
        residuals Res;
        hi Lev.
      let Press[i] = ssq(YY-((Trans-(Res/(1-Lev)))**(1/c)))
    ENDIF
ENDDO

#***************************************************************************
#   Plot log-likelihood versus lambda, PRESS versus lambda

let lo = round(lo,2)
let hi = round(hi,2)
let lam = round(lam,2)
copy lo Junk
text Junk Junk
copy Junk lo
copy hi Junk
text Junk Junk
copy Junk hi
copy lam Junk
text Junk Junk
copy Junk lam

kkset g " , "
kkset h "("
kkset j ")"
kkset k "Approximate 95% CI for Lambda: "
kkset l "Estimated Lambda: "

kkcat h lo h
kkcat h g h
kkcat h hi h
kkcat h j h

kkcat k h k
kkcat l lam l

copy lo Junk
numeric Junk Junk
copy Junk lo

copy hi Junk
numeric Junk Junk
copy Junk hi

copy lam Junk
numeric Junk Junk
copy Junk lam

copy l T1
copy k T2

brief
write T1
write T2
brief 0

title
layout;
title 'Box-Cox Power Transformation Analysis';
wtitle 'Box-Cox Power Transformation Analysis'.

plot LL*lambda;
 connect;
  color 4;
 refe 1 lo hi;
  type 2;
  color 2;
 grid 1;
  color 14;
 grid 2;
  color 14;
 marker lam e;
  color 2;
 axlabel 2 "Log-Likelihood";
 nodtitle;
 footnote l;
 footnote k;
 figure 0 0.5 0 1.

plot Press*Mesh;
 connect;
  color 4;
 grid 1;
  color 14;
 grid 2;
  color 14;
 axlabel 2 "PRESS";
 axlabel 1 "Lambda";
 nodtitle;
 figure 0.5 1 0 1.

endlayout
notitle

#***************************************************************************
#   Cook-Wang Index Plot

IF influence
  let m = (sum(Y1))/(count(YY))
  count YY a
  set Col
  a(1)
  end
  diagonal Col m2
  copy Col XX.1-XX.n m3
  transpose m3 m4
  multiply m3 m1 m5
  multiply m5 m4 m5
  subtract m5 m2 m5
  diagonal m5 Col1
  let Col1 = Col1**(-0.5)

    IF lam = 0
      let Col2 = (expo((sum(Y1))/(count(YY))))*Y1
      let Col3 = (Col2)*((Y1/2)-m)
    ELSE
      let Col2 = ((YY**lam)-1)/(lam*(expo((lam-1)*m)))
      let Col3 = (((YY**lam)*((lam*Y1)-(lam*m)-1))+(lam*m)+1)/((lam**2)*(expo((lam-1)*m)))
    ENDIF

  copy Col2 m6
  copy Col3 m7
  multiply m5 m7 m8
  multiply m5 m6 m9
  copy m8 Col4
  copy m9 Col5
  let Col4 = (Col4*Col1)/(((ssq(Col4))/(count(YY)-n-1))**0.5)
  let Col5 = (Col5*Col1)/(((ssq(Col5))/(count(YY)-n-1))**0.5)
  let Final = count(YY)*(loge(1+(((Col5*Col4)/((count(YY)-n-1)-(Col4*Col4)))**2)))

  set Case
  1:a
  end

title
tsplot Final;
  connect;
   color 4;
  symbol;
   color 2;
  axlabel 2 "Approximate Likelihood Distance";
  axlabel 1 "Case Number";
  title 'Cook-Wang Index Plot';
  grid 1;
   color 14;
  grid 2;
   color 14;
  wtitle 'Cook-Wang Index Plot'.
notitle

ENDIF

#***************************************************************************
#   Storage

IF bcstore
  copy LL Lambda St1 St2
ENDIF

IF infstore
  copy Final Case St3 St4
ENDIF

IF presstore
  copy Press Mesh St5 St6
ENDIF

#***************************************************************************

brief 2

title
Endmacro