Yet Another Blog in Statistical Computing

I can calculate the motion of heavenly bodies but not the madness of people. -Isaac Newton

Modeling Heteroscedasticity Directly in NLS

****** nonlinear least square regression without heteroscedasticity ******;
proc nlmixed data = data.sme tech = congra;
  where y > 0  and y < 1;
  parms b0 =  1.78 b1 = -0.01 b2 = -0.43 b3 = -0.11 b4 = -2.93
        b5 =  0.01 b6 = -0.01 b7 = -0.17; 
  xb = b0 + b1 * x1 + b2 * x2 + b3 * x3 + b4 * x4 +
       b5 * x5 + b6 * x6 + b7 * x7;
  mu = 1 / (1 + exp(-xb));
  lh = pdf('normal', y, mu, s);
  ll = log(lh);
  model y ~ general(ll);
run;
/*
             Fit Statistics
-2 Log Likelihood                 -264.5
AIC (smaller is better)           -246.5
AICC (smaller is better)          -246.3
BIC (smaller is better)           -201.4

                       Standard
Parameter   Estimate      Error     DF   t Value   Pr > |t|    
b0            1.7813     0.3400   1116      5.24     <.0001   
b1          -0.01203    0.02735   1116     -0.44     0.6602  
b2           -0.4305     0.1437   1116     -3.00     0.0028    
b3           -0.1147    0.02287   1116     -5.02     <.0001   
b4           -2.9302     0.4657   1116     -6.29     <.0001    
b5          0.004095   0.001074   1116      3.81     0.0001   
b6          -0.00839   0.002110   1116     -3.98     <.0001    
b7           -0.1710     0.1977   1116     -0.87     0.3871   
s             0.2149   0.004549   1116     47.24     <.0001  
*/
 
****** nonlinear least square regression with heteroscedasticity ******;
proc nlmixed data = data.sme tech = congra;
  where y > 0 and y < 1;
  parms b0 =  1.78 b1 = -0.01 b2 = -0.43 b3 = -0.11 b4 = -2.93
        b5 =  0.01 b6 = -0.01 b7 = -0.17 s  =  0.21
        a1 = -0.13 a2 = -15.62 a3 = 0.09 a4 = -1.27 a5 = 0.01
        a6 = -0.02 a7 =  0.47;
  xb = b0 + b1 * x1 + b2 * x2 + b3 * x3 + b4 * x4 +
       b5 * x5 + b6 * x6 + b7 * x7;
  xa = a1 * x1 + a2 * x2 + a3 * x3 + a4 * x4 +
       a5 * x5 + a6 * x6 + a7 * x7;
  mu = 1 / (1 + exp(-xb));
  si = (s ** 2 * (1 + exp(xa))) ** 0.5;       
  lh = pdf('normal', y, mu, si);
  ll = log(lh);
  model y ~ general(ll);
run;
/*
             Fit Statistics
-2 Log Likelihood                 -325.9
AIC (smaller is better)           -293.9
AICC (smaller is better)          -293.4
BIC (smaller is better)           -213.6

                       Standard
Parameter   Estimate      Error     DF   t Value   Pr > |t|    
b0            2.0343     0.3336   1116      6.10     <.0001     
b1          0.003764    0.02408   1116      0.16     0.8758    
b2          -0.08544     0.1501   1116     -0.57     0.5693     
b3           -0.1495    0.02263   1116     -6.61     <.0001    
b4           -2.6251     0.4379   1116     -6.00     <.0001    
b5          0.003331   0.001115   1116      2.99     0.0029    
b6          -0.00644   0.001989   1116     -3.24     0.0012     
b7           -0.1836     0.1938   1116     -0.95     0.3436    
s             0.1944   0.005067   1116     38.35     <.0001    
a1           -0.1266     0.3389   1116     -0.37     0.7088     
a2          -15.6169     5.2424   1116     -2.98     0.0030     
a3           0.09074    0.03282   1116      2.76     0.0058     
a4           -1.2681     3.8044   1116     -0.33     0.7389    
a5          0.007411   0.005267   1116      1.41     0.1597    
a6          -0.01738    0.01527   1116     -1.14     0.2550    
a7            0.4805     1.1232   1116      0.43     0.6689    
*/
Advertisements

Written by statcompute

October 7, 2012 at 12:46 am

%d bloggers like this: