R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("multinomRob-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('multinomRob')
Loading required package: rgenoud
Loading required package: MASS
Loading required package: mvtnorm
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "multinomRob"
> 
> ### * multinomRob
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: Robust Multinomial Regression
> ### Title: Multinomial Robust Estimation
> ### Aliases: multinomRob
> ### Keywords: robust models regression
> 
> ### ** Examples
> 
> # make some multinomial data
> x1 <- rnorm(50);
> x2 <- rnorm(50);
> p1 <- exp(x1)/(1+exp(x1)+exp(x2));
> p2 <- exp(x2)/(1+exp(x1)+exp(x2));
> p3 <- 1 - (p1 + p2);
> y <- matrix(0, 50, 3);
> for (i in 1:50) {
+   y[i,] <- rmultinomial(1000, c(p1[i], p2[i], p3[i]));
+ }
> 
> # perturb the first 5 observations
> y[1:5,c(1,2,3)] <- y[1:5,c(3,1,2)];
> y1 <- y[,1];
> y2 <- y[,2];
> y3 <- y[,3];
> 
> # put data into a dataframe
> dtf <- data.frame(x1, x2, y1, y2, y3);
> 
> ## Set parameters for Genoud
> zz.genoud.parms <- list( pop.size             = 1000,
+                         wait.generations      = 10,
+                         max.generations       = 100,
+                         scale.domains         = 5,
+                         print.level = 0
+                         )
> 
> # estimate a model, with "y3" being the reference category
> # true coefficient values are:  (Intercept) = 0, x = 1
> # impose an equality constraint
> # equality constraint:  coefficients of x1 and x2 are equal
> mulrobE <- multinomRob(list(y1 ~ x1, y2 ~ x2, y3 ~ 0),
+                       dtf,
+                       equality = list(list(y1 ~ x1 + 0, y2 ~ x2 + 0)),
+                       genoud.parms = zz.genoud.parms,
+                       print.level = 3, iter=TRUE);

Equality constraints among parameters (after consolidation):
Equality constrained set 1 
outcome y1 regressor x1 
outcome y2 regressor x2 

Your Model (xvec):
                           y1 y2 y3
(Intercept)/(Intercept)/NA  1  1  0
x1/x2/NA                    2  2  0


multinomRob(): Grouped MNL Estimation
[1] "multinomMLE: hessian determinant: 1001661941978.6"
[1] "multinomMLE: OPG determinant: 197763281862155488"
MNL LQD Fit: 4.105606 
MNL Estimates:
                                   y1         y2 y3
(Intercept)/(Intercept)/NA 0.01084097 0.02569787  0
x1/x2/NA                   0.74739947 0.74739947  0

MNL SEs:
                                   y1        y2 y3
(Intercept)/(Intercept)/NA 0.03972336 0.1161921  0
x1/x2/NA                   0.14261554 0.1426155  0



multinomRob(): Calculating multinomial-t starting values.
Multinom-T LQD Fit (step  0 ): 1.285969 
Multinom-T Beta Estimates (step  0 ):
                                    y1           y2 y3
(Intercept)/(Intercept)/NA 0.002128599 -0.007719271  0
x1/x2/NA                   0.992721278  0.992721278  0
Multinom-T Beta SEs (step  0 ):
                                    y1          y2 y3
(Intercept)/(Intercept)/NA 0.013171167 0.010447475  0
x1/x2/NA                   0.009890758 0.009890758  0
Multinom-T Omega Estimates (step  0 ):
             y1           y2
y1 0.0043516170 0.0009569967
y2 0.0009569967 0.0032056389
Multinom-T DF (step  0 ): 1.145307 

multinomRob(): Using multinomial-t estimates as starting values.
Multinom-T LQD Fit: 1.285969 
Multinom-T Estimates:
                                    y1           y2 y3
(Intercept)/(Intercept)/NA 0.002128599 -0.007719271  0
x1/x2/NA                   0.992721278  0.992721278  0
Multinom-T DF: 1.145307 

multinomRob(): Starting Values 
                                    y1           y2 y3
(Intercept)/(Intercept)/NA 0.002128599 -0.007719271  0
x1/x2/NA                   0.992721278  0.992721278  0
multinomRob(): starting fit = 1.285969 


LQD Results:
                                     y1          y2 y3
(Intercept)/(Intercept)/NA 2.649726e-05 -0.01027215  0
x1/x2/NA                   1.001777e+00  1.00177737  0

LQD sigma: 1.246868 


(multinomTanh):
[1] "mGNtanh: number of Newton iterations 3"
[1] "mGNtanh: number of Newton iterations 2"
[1] "mGNtanh: number of Newton iterations 2"
[1] "mGNtanh: number of Newton iterations 2"
[1] "mGNtanh: number of Newton iterations 1"
[1] "mGNtanh: hessian determinant: 522418582130.271"
[1] "mGNtanh: OPG determinant: 484383404062.361"
[1] "mGNtanh: tanh sigma^2: 1.01233347617417"
Tanh Estimates
                                    y1           y2 y3
(Intercept)/(Intercept)/NA 0.003956141 -0.004350162  0
x1/x2/NA                   1.006994224  1.006994224  0

Tanh Sandwich SEs
                                   y1         y2 y3
(Intercept)/(Intercept)/NA 0.01346586 0.01102388  0
x1/x2/NA                   0.01033832 0.01033832  0

TANH sigma: 1.006148 


**********************************

Iteration: 1 
using best non-genoud LQD sigma. Provided by multinomial-t.

(iter 1 ): LQD Results:
                                    y1           y2 y3
(Intercept)/(Intercept)/NA 0.003958694 -0.004011233  0
x1/x2/NA                   0.994886340  0.994886340  0

(iter 1 ): LQD sigma: 1.251319 

WARNING: New fit is worse.  Stopping.
It may be a good idea to restart with a larger GENOUD population.
> summary(mulrobE, weights=TRUE);

Choice 1 : y1 Estimates and SE:
                Est SE.Sand t.val.Sand
(Intercept) 0.00396  0.0135      0.294
x1          1.01000  0.0103     97.400


Choice 2 : y2 Estimates and SE:
                 Est SE.Sand t.val.Sand
(Intercept) -0.00435  0.0110     -0.395
x2           1.01000  0.0103     97.400


Choice 3 : y3 Estimates and SE:
            Est SE.Sand t.val.Sand
NA            0       0        NaN
NA            0       0        NaN


LQD sigma: 1.246868 
TANH sigma: 1.006148 

Number of Observations: 50 
Number of observations with at least one zero weight: 5 
Number of zero weights: 9 
TANH: weights
   name weights:y1 weights:y2
1     1   0.000000          0
2     2   0.000000          0
3     3   0.000000          0
4     4   0.000000          0
5     5   0.511593          0
6     6   1.000000          1
7     7   1.000000          1
8     8   1.000000          1
9     9   1.000000          1
10   10   1.000000          1
11   11   1.000000          1
12   12   1.000000          1
13   13   1.000000          1
14   14   1.000000          1
15   15   1.000000          1
16   16   1.000000          1
17   17   1.000000          1
18   18   1.000000          1
19   19   1.000000          1
20   20   1.000000          1
21   21   0.922289          1
22   22   1.000000          1
23   23   1.000000          1
24   24   1.000000          1
25   25   1.000000          1
26   26   1.000000          1
27   27   1.000000          1
28   28   1.000000          1
29   29   1.000000          1
30   30   1.000000          1
31   31   1.000000          1
32   32   1.000000          1
33   33   1.000000          1
34   34   0.871220          1
35   35   1.000000          1
36   36   1.000000          1
37   37   1.000000          1
38   38   1.000000          1
39   39   1.000000          1
40   40   1.000000          1
41   41   1.000000          1
42   42   1.000000          1
43   43   1.000000          1
44   44   1.000000          1
45   45   1.000000          1
46   46   1.000000          1
47   47   1.000000          1
48   48   1.000000          1
49   49   1.000000          1
50   50   1.000000          1
> 
> 
> 
> cleanEx(); ..nameEx <- "rmultinomial"
> 
> ### * rmultinomial
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: rmultinomial
> ### Title: Random Number Generator for the Multinomial Distribution
> ### Aliases: rmultinomial
> ### Keywords: distribution multivariate
> 
> ### ** Examples
> 
>  rmultinomial(10, c(.3, .3, .4));
[1] 2 3 5
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  64.48 5.28 70.16 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')