N <- 100

myround <- function(x, dp)
{
  format <- paste("%.", dp, "f", sep="")
  sprintf(format, x)
}

# Make standardized variables with our desired correlations - this took a bit of trial and error
set.seed(1)
x <- scale(rnorm(N))
y <- scale(rnorm(N) + (x * .2383))
m <- scale(rnorm(N) + (x * .552) + (y * .311))

ryx <- cor(y, x)
cat("r(y,x)=", myround(ryx, 2), "\n", sep="")
rym <- cor(y, m)
cat("r(y,m)=", myround(rym, 2), "\n", sep="")
rxm <- cor(x, m)
cat("r(x,m)=", myround(rxm, 2), "\n", sep="")

reg <- lm(y ~ x + m)
betas <- reg$coefficients
tstats <- summary(reg)[["coefficients"]][, "t value"]

# First calculation of R-squared
# http://faculty.cas.usf.edu/mbrannick/regression/Part3/Reg2.html
Rsq <- (betas["x"] * ryx) + (betas["m"] * rym)

# Calculate standardized coefficients using only the correlation matrix.
# https://stats.stackexchange.com/questions/107597/is-there-a-way-to-use-the-covariance-matrix-to-find-coefficients-for-multiple-re
m.ryx <- .24
m.rym <- .32
m.rxm <- .52
# Manual R-squared calculation from http://faculty.cas.usf.edu/mbrannick/regression/Part3/Reg2.html
m.Rsq <- ((m.ryx ^ 2) + (m.rym ^ 2) - (2 * m.ryx * m.rym * m.rxm)) / (1 - (m.rxm ^ 2))
m.betas <- solve(matrix(c(1, m.rxm, m.rxm, 1), ncol=2), c(m.ryx, m.rym))

# https://www3.nd.edu/~rwilliam/stats1/x91.pdf
# Note: This code only works for two predictors.
# If we add more predictors, the calculation of K in the next line will not be sufficient.
K <- length(m.betas) - 1
se <- sqrt((1 - m.Rsq) / ((1 - (m.rxm ^ 2)) * (N - K - 1)))
tstats2 <- list("Intercept" = NA, "x" = m.betas[1] / se, "m" = m.betas[2] / se)

# Showing off now: We can calculate the ANOVA, then go from F ratios to t ratios.
# http://faculty.cas.usf.edu/mbrannick/regression/Part3/Reg2.html
## Fx <- (m.Rsq - m.rym^2) / ((1 - m.Rsq) / (N - K - 1))   #rym^2 is correct (R-squared without x).
## Fm <- (m.Rsq - m.ryx^2) / ((1 - m.Rsq) / (N - K - 1))
## tstats3 <- list("Intercept" = NA, "x" = sqrt(Fx), "m" = sqrt(Fm))

cat("Regressing Y on X and M\n")
cat("V1: numbers obtained from regressing the data points\n")
ivnames <- c("x", "m")
for (v1 in ivnames) {
  beta1 <- betas[v1]
  tv1 <- tstats[v1]
  p1 <- 2 * pt(-abs(tv1), N - 3)
  cat(v1, ": beta=", myround(beta1, 3), " t=", myround(tv1, 3), " p=", myround(p1, 3), "\n", sep="")
}

cat("V2: numbers obtained from the correlation matrix\n")
for (i in 1:length(ivnames)) {
  v2 <- ivnames[i]
  beta2 <- m.betas[i]
  tv2 <- tstats2[v2][[1]][[1]]
  p2 <- 2 * pt(-abs(tv2), N - 3)
  cat(v2, ": beta=", myround(beta2, 3), " t=", myround(tv2, 3), " p=", myround(p2, 3), "\n", sep="")
}