lm.perm<-function(y,x1,nperm=5000){
  obs.f<-anova(lm(y~x1))[, 4]
  bf<-matrix(NA,nrow=nperm,ncol=length(obs.f))
  #bf<-matrix(numeric(length(nperm)*length(obs.f)), nrow=nperm)
  bf[nperm,]<-obs.f
  
  for (i in 1:(nperm-1)){
    bf[i,]<-anova(lm(sample(y)~x1))[,4]
  }
  
  p.val<-apply(bf,2,fun<-function(x) sum(x>=x[length(x)])/nperm)
  
  return(list(p.val=p.val))
  
}

lm.perm(data$comb1,data$treatment)

m1=lm(data$comb1~data$treatment)
anova(m1)