Como calcular a variância do estimador de mínimos quadrados usando a decomposição QR em R?

Estou tentando aprender a decomposição do QR, mas não consigo descobrir como obter a variação do beta_hat sem recorrer aos cálculos matriciais tradicionais. Estou praticando com oiris conjunto de dados, e aqui está o que tenho até agora:

y<-(iris$Sepal.Length)
x<-(iris$Sepal.Width)
X<-cbind(1,x)
n<-nrow(X)
p<-ncol(X)
qr.X<-qr(X)
b<-(t(qr.Q(qr.X)) %*% y)[1:p]
R<-qr.R(qr.X)
beta<-as.vector(backsolve(R,b))
res<-as.vector(y-X %*% beta)

Obrigado pela ajuda!