Hello List and thanks in advance for all of your help,


I am trying implement a permutation test of a multinomial logistic
regression ('multinom' within the nnet package).  In the end I want to
compare the parameter estimate from my data to the distribution of
randomized parameter estimates.

I have figured out how to permute my dependent variable (MNNUM) x number of
times, apply multinomial logistic regression, to each permutation, and save
the results in a list.  Where I am stuck is figuring out how to take the
mean and SD of the coefficients from my list of regressions.  I know that
the coefficients are stored in the $wts slot of the model.

Below is what I have so far.  I am sure there are nicer ways to do this and
if you feel so inclined please suggest them.



#this is a function to permute the MNNUM column once
rand<- function(DF){
new.DF<-DF
new.DF$MNNUM<-sample(new.DF$MNNUM)
new.DF

}

#this function does one model I am interested in.

modeltree<-function(DF){
MLM.plot <- multinom(MN_fact  ~ Canpy  + mean_dbh  + num_beechoak  +
num_class5  + prop_hard , data=hfdata, trace=FALSE)
MLM.plot
}

# this replicates the 'rand' function and applies a model

resamp.funct<-function(DF,funct, n){
    list<-replicate(n,rand(DF), simplify = FALSE)
    sapply(list, funct, simplify = FALSE)
}


#So if I paste below:

l<-resamp.funct(hfdata, modeltree, 3)

# I get


> l<-resamp.funct(hfdata, modltree, 3)

> l
[[1]]
Call:
multinom(formula = MN_fact ~ Canpy + mean_dbh + num_beechoak +
    num_class5 + prop_hard, data = hfdata, trace = FALSE)

Coefficients:
         (Intercept)       Canpy    mean_dbh num_beechoak num_class5
prop_hard
none     -11.1845028 0.063880939  0.08440340   -0.7050239 -0.0998379
6.894522
sabrinus -10.6848488 0.055157318  0.19276777   -0.6441996  0.1219245
3.325704
volans    -0.2481854 0.004410597 -0.02710102   -0.1061700 -0.1858376
2.495856

Residual Deviance: 163.7211
AIC: 199.7211

[[2]]
Call:
multinom(formula = MN_fact ~ Canpy + mean_dbh + num_beechoak +
    num_class5 + prop_hard, data = hfdata, trace = FALSE)

Coefficients:
         (Intercept)       Canpy    mean_dbh num_beechoak num_class5
prop_hard
none     -11.1845028 0.063880939  0.08440340   -0.7050239 -0.0998379
6.894522
sabrinus -10.6848488 0.055157318  0.19276777   -0.6441996  0.1219245
3.325704
volans    -0.2481854 0.004410597 -0.02710102   -0.1061700 -0.1858376
2.495856

Residual Deviance: 163.7211
AIC: 199.7211

[[3]]
Call:
multinom(formula = MN_fact ~ Canpy + mean_dbh + num_beechoak +
    num_class5 + prop_hard, data = hfdata, trace = FALSE)

Coefficients:
         (Intercept)       Canpy    mean_dbh num_beechoak num_class5
prop_hard
none     -11.1845028 0.063880939  0.08440340   -0.7050239 -0.0998379
6.894522
sabrinus -10.6848488 0.055157318  0.19276777   -0.6441996  0.1219245
3.325704
volans    -0.2481854 0.004410597 -0.02710102   -0.1061700 -0.1858376
2.495856

Residual Deviance: 163.7211
AIC: 199.7211

	[[alternative HTML version deleted]]

