Bootstrapping linear regression models in R (wild and pairs)

Question

I am trying to show the differences between the residual, wild and pairs bootstrap on a regression model in R. I have some small issues with my code for the pairs bootstrap and I would like to use the wild bootstrap on the same data to get some estimates but I am unsure of how to do this, is this possible?

so far I have successfully coded this for the residual (using the ToothGrowth data in R):

library(nptest)
xmat<- cbind(1, ToothGrowth$dose)
xinv<- solve(crossprod(xmat)) %*% t(xmat)
fit<- xmat %*% xinv %*% ToothGrowth$len
data<- list(fit=fit, resid=ToothGrowth$len - fit, xinv = xinv, x=ToothGrowth$dose)
statfun<- function(x, data) {
  ynew<- data$fit+data$resid[x]
  coef<-as.numeric(data$xinv %*% ynew)
  names(coef)<-c("Intercept", "len")
  coef
}

jackfun<-function(x, data){
  ynew<-data$fit[x]+data$resid[x]
  xmat<-cbind(1, data$x[x])
  xinv<-solve(crossprod(xmat)) %*% t(xmat)
  coef<- as.numeric(xinv %*% ynew)
  names(coef)<- c("Intercept", "len")
  coef
}

npbs<-np.boot(x=1:nrow(ToothGrowth), statistic= statfun,
                       data=data, jackknife = jackfun)
npbs        

Which returns nice estimates, standard error, bias and CIs. Then for the pairs bootstrap I have done this:

library(L1pack)
library(boot)
fit.model<-function(ToothGrowth)
  {fit<-glm(len~dose, data = ToothGrowth)
  l1<-l1fit(ToothGrowth$len, ToothGrowth$dose)
  c(coef(fit),coef(l1)) }
ToothGrowth.fun<-function(ToothGrowth, i) fit.model(ToothGrowth[i,])
tooth.boot1<-boot(ToothGrowth, ToothGrowth.fun, R=1000)
tooth.boot1

Which does give nice estimates, however 1) it has an error about non-unique solution possible which I would like to resolve. It also returns 4 estimates whereas I am only interested in the first two, but when I remove the function part, it does not run.

I am unsure on how to start the code for the wild bootstrap.