qtest.src

没有说明

/*
** qtest.src - Wald test of Logit and probit model
** (C) Copyright 1988-1998 by Aptech Systems, Inc.
** All Rights Reserved.
**
** This Software Product is PROPRIETARY SOURCE CODE OF APTECH
** SYSTEMS, INC.    This File Header must accompany all files using
** any portion, in whole or in part, of this Source Code.   In
** addition, the right to create such files is strictly limited by
** Section 2.A. of the GAUSS Applications License Agreement
** accompanying this Software Product.
**
** If you wish to distribute any portion of the proprietary Source
** Code, in whole or in part, you must first obtain written
** permission from Aptech Systems.
**
**> qtest
**
**
**  Purpose:  To test linear hypotheses for the quantal response models.
**            The details of the model are presented in the remarks section.
**
**  Format:   { Wald } = QTEST(indvars,b,vc,teststr);
**
**  Input:
**            indvars    Kx1 character vector, names of independent variables.
*                          - or -
**                       Kx1 numeric vector, indices of independent variables.
**
**                       The program adds one variable for the constant term.
**
**
**            b    an NPARM=(NCAT-1)*(K+1) vector of parameter estimates in
**                 the order:  intercepts|var1|var2|...varK.  For each
**                 variable the parameters are in the order comparing the
**                 first category to NCAT, the second to NCAT, ... to
**                 NCAT-1 to NCAT.  See below for details.
**
**            vc    NPARMxNPARM variance covariance matrix for the parameters
**                  in b.
**
**        teststr    String of the test statement.  For example, if you want
**                   to test whether the coefficient of variable x1 is equal
**                   to x2's,  you should write:
**
**              teststr = "x1 = x2";
**
**  Output:
**       Wald   Wald statistics of linear hypotheses test.
**
**  Globals:
**
**            __output    global scalar, default 1.
**
**                        if 1, sends results to the output device (including
**                        the screen).
**
**                        if 0, no information is sent to output.
**
**  Remarks:
**      You must run the model to be tested first to get input for this
**      test program, then you can do any number of tests you want.
**
**
**  Example:
**
**      library quantal;
**      #include quantal.ext
**      quantset;
**
**      output file = test.out reset;
**      __miss = 1;
**      dsn = "aldnel";
**      dv = { a };
**      iv = { gpa, tuce, psi };
**      test1 = " gpa + 2 tuce = 0";
**      test2 = " tuce = psi";
**      { vnames,b,vc,n,pct,mn,sd,fit,df,tol } = logit(dsn,dv,iv);
**      { Wald1 } = QTEST(vnames,b,vc,test1);
**      { Wald2 } = QTEST(vnames,b,vc,test2);
**      output off;
**
**  Library:  quantal
**
**  See Also: logit, ordered, probit, psnreg
*/

#include gauss.ext
#include quantal.ext

proc(1) = Qtest(vnames,b,vc,strng);

  local Wald,CHIFIG,err,z,R,str,invar,nct,rhnum,
        errmsg1,nvar,rowvc,i,nmcma,rmtrx,temp,ex;

  if __output;
     print;
     print "=================================================="\
       "===========================";
     print;
     print"The Hypothesis to be tested is: ";
     print $strng;
     print;
     print;
     print "*****************  Results from Linear Hypothesis Testing  "\
           "*****************";
     print;
  endif;

  clear nmcma;
  temp = vals(strng);
  strng = chrs(packr(miss(miss(miss(temp,10),13),32)));
  temp = vals(strng);
  ex = (temp ./= 44);
  ex = packr(missex(temp,ex));
  if not ismiss(ex);
      nmcma = rows(packr(ex));
  endif;
  clear temp;
  invar = trimr(vnames, 1, 0);
  nvar = rows(invar);
  rowvc = rows(vc);
  nct = rowvc/nvar;
  if round(nct) /= nct;
      errmsg1 = "ERROR: Incorrect arguments passed to this procedure";
      goto errout(error(7));
  endif;
  rhnum = nvar*ones(nct,1);
  invar = reshape(invar,rowvc,1);
  if nct == 1;
      { r,z } = rmatrix(strng,invar);
      if scalerr(r);
          goto errout(r);
      endif;
  else;
      { r, z } = srmatrix(strng,invar,rhnum);
      if scalerr(r);
         goto errout(r);
      endif;
      if nmcma > 0;
           i = 1;
           do until i > nmcma + 1;
               rmtrx = reshape(r[i,.],nct,nvar);
               rmtrx = reshape(rmtrx',rowvc,1);
               r[i,.] = rmtrx';
               i = i + 1;
          endo;
      endif;
  endif;
  if scalmiss(b);
      errmsg1 = "ERROR:  The estimated coefficients passed to"\
            " this procedure is a missing value.";
      goto errout(error(10));
  endif;

  local scale, rvci, s;
  scale = sqrt(diag(vc));   /*  scale the matrices */
  vc = vc./scale./scale';
  R = R.*scale';
  b = b./scale;
  rvci = invswp(R*vc*R');
  s = sumc(diag(rvci) > 0);
  if s > 0;
      Wald = (R*b-z)'*rvci*(R*b-z);
  else;
      errmsg1 = "ERROR: Wald test failed";
      goto errout(0);
  endif;

/*
  Remark: in Shazam, they set z = 0; even if user supplies z /= 0
  so the result is different if z /= 0
*/

  CHIFIG = cdfchic(Wald,rows(R));

  if __output;
     str = ftos(s,"       Wald Chi-SQ(%*.*lf) statistic = ",1,0);
     str = strsect(str,1,34) $+ ftos(Wald,"%*.*lf",7,3);
     print str;;
     print ftos(CHIFIG,"      Prob. = %*.*lf",6,3);
     print;
     if s < rows(R);
         print "Warning: rank of statistic less than order of cov"\
           "ariance matrix";
     print;
     endif;
     print "*****************************************************"\
           "*************************";
     print;
  endif;

  retp(Wald);

errout:
    pop err;
    if not trapchk(1);
       errorlog errmsg1;
       print;
       end;
    endif;
    retp(err);

endp;