I believe that this code will work (...for very small) samples, but let some
correct me if there is something wrong.

require(logistf);require(combinat)
permY<-permn(data$y)
ntimes<-length(permY)
results<-matrix(nrows=ntimes,ncols=number_of_coefficients)
for(i in 1:ntimes){
results[i,]<-logistf(unlist(permY[i])~factor_A+factor_B,data=data)}
observed<-logistf(y~factor_A+factor_B,data=data)$coefficients
# Exact p-values will be:
sum(results[,1]>=observed[1])/ntimes # One for intercept
sum(results[,2]>=observed[2])/ntimes
sum(results[,3]>=observed[3])/ntimes
-- 
Miłego dnia

	[[alternative HTML version deleted]]

