Subject: CCL:constrained multiple linear regression

From: [email-masked]
Date: Mon, 2 Dec 2002 19:08:14 -0600

CCL,

        I am looking for software (preferably free or inexpensive, code is 
OK) that will allow me to perform
constrained multiple linear regression.  For example, I have an equation 
of the form:

        Y = a*X1 + b*X2 + c*X3 + d

        where: Y is the dependent variable
                  X1 - X3 are independent variables (more than 1)
                a, b, c are coefficients
                        d is a constant

        I must be able to constrain the range of values for a, b, c, and d 
during the fitting.

        Please note that programs such as JMP, KalediaGraph, GraphPad 
Prism, and several other packages do 
NOT allow the user to do this, especially for more than one independent 
variable.

        Can someone suggest a program or code to do this?

        Thank you,
        Jim Metz

James T. Metz, Ph.D.
Research Investigator Chemist

GPRD R46Y AP10-2
Abbott Laboratories
100 Abbott Park Road
Abbott Park, IL  60064-6100
U.S.A.

Office (847) 936 - 0441
FAX    (847) 935 - 0548

[email-masked]

---

From: Jan Labanowski <[email-masked]>
Date: Mon, 2 Dec 2002 22:20:57 -0500 (EST)

You need a program which does regression through origin (i.e., no
intercept regression. Check references (sorry, it was ages ago when I looked
at it):
     K.W. Smillie - "An introduction to regression and correlation"
        Academic Press, NY 1966 
     S.R. Searle - "Linear models",  J.Wiley & Sons, NY 1971
   
I once write my own routine to do it, but it is in a terribly ancient Fortran
(yes... remember VAXes? in the early 80-ties?).
You have a regression
         Y = X*b + e   
When freeft = 1, then we fit an intercept and the X contains a column of 1's
When freeft = 0, then we have regression throught origin and X does not have
                 a column of 1's

(Note, you can always subtract a constant from all values of Y if you want to
fix intercept at some value. Likewise, if you want to fix some b, you
just subtract the b(j)*X(i,j) from Y(i) )
        Y - column vector of "experimental values"
        X - matrix storing values of independent variables:
            X(i,j) - value j-th independent variable for i-th experimental
        b - column vector of regression coefficients
            b(i) - regresion coefficients corresponding to i-th independent
                   variable
        e - column vector of residuals
        Nexp - number of experimental points
        Nparam - number of fitted regression coefficients

If you have statistical weights, W (i.e., reciprocals of variance for each
experimental point), calculate their mean
    
Wmean = (sum of Wi)/Nexp
and calculate normalized roots of weights as:

Wroot(i) = sqrt(W(i)/Wmean)

Then rescale Y and X by Wroots as:

Y(i) = Y(i) * Wroot(i)    (for i = 1, Nexp)

X(i,j) = X(i,j) * Wroot(i) (for i = 1, Nexp and for j = 1, Nparam)

Then calculate:

C             XTX - on output inverse of X'X matrix (inside the routine,
C                   before inversion, it stores elements of X'X)
C             XTY - X'Y matrix
C             YTY - Y'Y - sum of squares of Yexp(i)
C             rXiXj - array of correlation coefficients between independent
C                   variables Xi and Xj
C             rYXi - array of correlation coefficients between dependent var-
C                  iable Y and independent variable Xi
C             SSE - sum of squares of residuals
C                   SSE = sum of (Yexp(i) - Ytheor(i))**2  (i.e Y'Y - b'X'Y)
C                   both for freeft = 0  &  freeft = 1
C
C             SSR - sum of sqyares about regression
C                   SSR = sum of (Ytheor(i) - Ymean)**2
C                        or b'X'Y - Nexp*Ymean*Ymean  when freeft = 1
C                   SSR = b'X'Y  when freeft = 0
C
C             SST - total sum of squares
C                   SST = (Yexp(i) - Ymean)**2
C                       or Y'Y - Nexp*Ymean*Ymean  when freeft = 1
C                   SST = Y'Y   when  freeft = 0
C
C             MSE - mean square of errors  { MSE = SSE/(Nexp-Nparam)  }
C
C             MSR - mean square of regression
C                   MSR = SSR/Nparam for freeft = 0
C                   MSR = SSR/(Nparam - 1) for freeft = 1
C
C             b   - regression coefficients b = inv(X'X)*X'Y
C                   inv(X'X) is the inverse of X'X 
C
C             seb - standard errors of regression coefficients
C                   seb[i] = SQRT(MSE/inv(X'X[i,i]))
C                          where  inv(X'X[i,i]) is the i,i-element of the
C                                 inverse of X'X matrix 
C
C             seYth - standard errors of regression values
C                     seYth[i] = SQRT(MSE*Xi'*inv(X'X)*Xi) where
C                     Xi is the i-th column of X matrix and
C                     inv(X'X) is the inverse of X'X matrix
C
C             Rmul - multiple regression coefficient that is SQRT(Rmul2) or
C                    in case of Rmul2 < 0 (sometimes in no-intercept model)
C                    Rmul is set to 0
C
C             Rmul2 - Square of multiple correlation coefficient
C                     Rmul2 = SSR/SST
C
C             Xmean - array of means of independent variables
C
C             Ymean - mean of Yexp  (and Ytheor)
C
C             F -  F statistics         F = MSR/MSE
C
C             t - t statistics for each regression coefficient for testing
C                  hypothesis that bstar(i) = 0
C                     t(i) = b[i]/seb[i]
C
C The ugly Fortan follows

       INTEGER*2 Nexp, Nparam, df, i, j, k, l, errflg, mode,freeft
       REAL*8    X(MAXEXP,MAXPAR), Yexp(MAXEXP), Ytheor(MAXEXP),
     1           W(MAXEXP), Wroots(MAXEXP), Wmean, t(MAXPAR),
     2           rXiXj(MAXPAR,MAXPAR), rYXi(MAXPAR), Xmean(MAXPAR),
     3           XTX(MAXPAR,MAXPAR), XTY(MAXPAR), SSE, SSR, SST,
     4           seYth(MAXEXP), seb(MAXPAR), b(MAXPAR), Rmul2,
     5           Rmul, F, MSR, MSE, Ymean, YTY, y, s, s1, s2, aux

C Read in the Yexp, X, W and freeft (find the return statements and
C change them to STOP) or make a subroutine out of this code
C ....................
C

C
C    Calculate Wmean for relative weights
C
       Wmean = 0.0D0
       DO 10 i = 1, Nexp
          Wmean = Wmean + W(i)
10     CONTINUE

       Wmean = Wmean/Nexp
C
C    Rescale Yexp and X for weighted least squares
C

15     DO 20 i = 1, Nexp
          Wroots(i) = SQRT(W(i)/Wmean)
20     CONTINUE

       DO 40 i = 1, Nexp
          s = Wroots(i)
          Yexp(i) = Yexp(i) * s
          DO 30 j = 1, Nparam
             X(i,j) = X(i,j) * s
30        CONTINUE
40     CONTINUE

C
C    Calculate Y'Y and Ymean
C
       YTY = 0.0D0
       Ymean = 0.0D0
       DO 42 i = 1, Nexp
          y = Yexp(i)
          YTY = YTY + y*y
          Ymean = Ymean + y
42     CONTINUE
       Ymean = Ymean/Nexp

C
C    Calculate Xmean
C
       DO 44 i = 1, Nparam
          s = 0.0D0
          DO 43 j = 1, Nexp
             s = s + X(j,i)
43        CONTINUE
          Xmean(i) = s/Nexp
44     CONTINUE

C
C    Calculate X'X
C
       DO 70 i = 1, Nparam
          DO 60 j = i, Nparam
             s = 0.0D0
             DO 50 k = 1, Nexp
                s = s + X(k,i)*X(k,j)
50           CONTINUE
             XTX(i,j) = s
             IF(i .ne. j) XTX(j,i) = s
60        CONTINUE
70     CONTINUE

C
C    Calculate X'Y
C
       DO 90 i = 1, Nparam
          s = 0.0D0
          DO 80 j = 1, Nexp
             s = s + X(j,i) * Yexp(j)
80        CONTINUE
          XTY(i) = s
90     CONTINUE

C
C    Calculate correlation coefficients rXiXj and rYXi
C
       DO 92 i = 1, Nparam
          s1 = Nexp * Xmean(i)**2
          DO 91 j = i, Nparam
             s2 =  Nexp * Xmean(j)**2
             s = XTX(i,j) - Nexp * Xmean(i) * Xmean(j)
             IF (ABS(XTX(i,i) - s1) .LT. 0.00001D0*s1 .OR.
     1           ABS(XTX(j,j) - s2) .LT. 0.00001D0*s2) THEN
                 rXiXj(i,j) = 0.0D0
             ELSE
                 aux = (XTX(i,i)-s1)*(XTX(j,j)-s2)
                 IF(aux .le. 1.0d-30)THEN
                   errflg = 30
                   return
                 end if
                 rXiXj(i,j) = s / SQRT((XTX(i,i)-s1)*(XTX(j,j)-s2))
             ENDIF
             IF(i .NE. j) rXiXj(j,i) = rXiXj(i,j)
91        CONTINUE
92     CONTINUE

       DO 93 i = 1, Nparam
          s = XTY(i) - Nexp * Xmean(i) * Ymean
          s1 = Nexp * Xmean(i)**2
          s2 = Nexp * Ymean*Ymean
          IF (ABS(XTX(i,i) - s1) .LT. 0.00001D0*s1 .OR.
     1        ABS(YTY - s2) .LT. 0.00001D0*s2) THEN
              rYXi(i) = 0.0D0
          ELSE
              aux = (XTX(i,i)-s1)*(YTY-s2)
              if(aux .le. 1.0D-30)THEN
                errflg = 30
                return
              end if
              rYXi(i) = s / SQRT((XTX(i,i)-s1)*(YTY-s2))
          ENDIF
93     CONTINUE

C
C    Invert X'X
C
       CALL syminv(Nparam, XTX, s, errflg)
       IF(errflg .ne. 0)GOTO 999

C
C    Calculate regression coefficients b = inv(X'X)*X'Y
C
       DO 110 i = 1, Nparam
          s = 0.0D0
          DO 100 j = 1, Nparam
             s = s + XTX(i,j) * XTY(j)
100       CONTINUE
          b(i) = s
110    CONTINUE

C
C    Calculate Ytheor (values of y given by regression) Ytheor = X*b
C
       DO 130 i = 1, Nexp
          s = 0.0D0
          DO 120 j = 1, Nparam
             s = s + b(j) * X(i,j)
120       CONTINUE
          Ytheor(i) = s
130    CONTINUE

C
C    Calculate SST
C
       IF(freeft .eq. 1)THEN
          SST = YTY - Nexp*Ymean*Ymean
       ELSE
          SST = YTY
       END IF

C
C    Calculate sum of squares of residuals (SSE)
C
       SSE = YTY
       DO 140 i = 1, Nparam
          SSE = SSE - XTY(i)*b(i)
140    CONTINUE

       IF(SSE .lt. -5.0d-10*YTY)THEN
          errflg = 2
          RETURN
       ELSE IF(SSE .lt. 0.0D0)THEN
         SSE = 0.0D0
       END IF

C
C    Calculate Rmul2,  and SSR (sum of squares due
C    to regression )
C
       df = Nexp - Nparam

       SSR = SST - SSE

       if(abs(SST) .gt. 1.0D-30)then
         Rmul2 = SSR/SST
       else
         Rmul2 = 0.0D0
       end if

C    mean square due to regression (MSR) and mean square of error (MSE)

       IF(freeft .eq. 0)THEN
          MSR = SSR/Nparam
       ELSE
          MSR = SSR/(Nparam - 1)
       END IF

       MSE = SSE/df

C    Calculate F statistic (F), multiple regression coefficient (Rmul)

       if(ABS(MSE) .gt. 1.0D-30)then
         F = MSR/MSE
       else
         F = 0.0D0
       end if

       IF(Rmul2 .gt. 0.0D0)THEN
          Rmul = SQRT(Rmul2)
       ELSE
          Rmul = 0.0D0
       END IF

C    Calculate t-statistic and standard error for each coefficient

       DO 150 i = 1, Nparam
          IF(XTX(i,i)*MSE .gt. 1.0D-30)THEN
            s = SQRT(MSE * XTX(i,i))
            seb(i) = s
            t(i) = b(i)/s
          ELSE
            seb(i) = 0.0D0
            t(i) = 0.0D0
          END IF
150    CONTINUE

C    Calculate standard error for each value of Ytheor

       DO 180 i = 1, Nexp
          s1 = 0.0D0
          DO 170 j = 1, Nparam
             s = 0.0D0
             DO 160 k = 1, Nparam
                s = s + XTX(j,k) * X(i,k)
160          CONTINUE
             s1 = s1 + s * X(i,j)
170       CONTINUE
          if(MSE*s1 .le. 1.0D-30)then
            seYth(i) = 0.0D0
          else
            seYth(i) = SQRT(MSE*s1)
          end if
180    CONTINUE

C    Restore original values of Yexp and X

       DO 200 i = 1, Nexp
          s = 1.0D0/Wroots(i)
          Yexp(i) = Yexp(i) * s
          DO 190 j = 1, Nparam
             X(i,j) = X(i,j) * s
190        CONTINUE
200    CONTINUE

On Mon, 2 Dec 2002 [email-masked] wrote:

> Date: Mon, 2 Dec 2002 19:08:14 -0600
> From: [email-masked]
> To: [email-masked]
> Cc: [email-masked]
> Subject: CCL:constrained multiple linear regression
> 
> CCL,
> 

Jan K. Labanowski            |    phone: 614-292-9279,  FAX: 614-292-7168
Ohio Supercomputer Center    |    E-mail: [email-masked] 
1224 Kinnear Rd,             |    http://www.ccl.net/~jkl
Columbus, OH 43212-1163      |    http://www.ccl.net/    http://asdn.net/