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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
