use R plz to write it

Linear models have the form Y = X beta + epsilon, where Y is the vector of length n, beta is the vector containing the p parameters, X is the design matrix (i. e. the matrix with the predictors as columns) of dimension n times p and epsilon represents Gaussian noise with mean zero and variance sigma^2middotThe parameter estimates are beta^cap = (X' X)^-1X' Y, (X'X stands for the matrix multiplication of the transpose of X with X, and (X' X)^-1 is the inverse of that product), an estimate for the variance is sigma^cap2 = (Y - xbeta^cap)(Y - Xbeta^cap)/(n - p), an estimate for the covariance matrix of the parameter estimates is cov^hat(beta^hat) = sigma^hat2, (X'X)^-1, the projection ("hat") matrix is P = X (X'X)^-1X' the fitted values are Y^hat = PY, and the residuals are Y^hat - Y. Write a function that takes the response vector Y and the matrix of covariates X as input, and returns a list of the following: beta, the vector of least squares estimates, sigma, the residual standard error, var. beta, the covariance matrix of the least squares estimates, fitted, the vector of fitted values, residuals, the vector of residuals. Further, put in an option to return hat, the projection matrix, upon request. The default should be to not return it. To fit an intercept, the elements in the first column of X have to be equal to one, so your function should also have an option to add a vector of ones to the matrix with the predictors. Your function should check whether or not X'X is invertible, and stop if it is not. Further, give warning when X'X is close to singular