How can I create meaningful / interesing data following multiple linear models?

Question

I would like to create many(!) examples for exercices, where I would like to control various parameters such as the coefficients, their uncertainty, t-values and p-values of the linear model.
I know, there is a lot of data in the web, but it is really difficult to find appropriate data to a specific question / example and if you need a new example you need to search the web again and again.
So I thought, I'll create data on my own - but it seems, this is really difficult. What I do not want

• Decrease n until t-/ p-values improve ("leading to no data")
• Increase sigma in rnorm until t-/ p-values improve ("an overall increase in sigma of a factor of 10 would do the job, but then no linear model is left." See the example below.)

Below you can find what I have using R. Is there a way to "improve" the result? In my concrete case I would like data which can be checked by plotting the data (about 100 to 1000 points, R's Std Error small) and p-values ranging from say 1e-5 to 0.8.

Edit: Thanks to Mickybo Yakari's answer, the situation improved a lot ($x_i$ values are sampled according to a multivariate gaussion distribution), but it would be great, if I could also "control" the p-values. The example below shows e.g. a Pr(>|t|)value for (Intercept) of 0.00016 but I would like this parameter to be more significant.

Is there any way to obtain, what I want?

#' Generate sparse precision matrix (Mickybo Yakari's answer fixed correlations)
#'
#' @param dimension An integer, the number of rows of the precision matrix.
#' @param upper A numeric in (0,1) specifying the range of allowed non-zero entries.
#' @param seed An integer, the random seed.
#'
#' @return A precision matrix
generate.sparse.precision.matrix <- function(dimension, upper, seed) {
  matrix <- matrix(rep(0,dimension*dimension), ncol = dimension)
  set.seed(seed)
  vec <- runif(n = dimension^2, min = 0, max = 1)

  for (i in 1:dimension) {
    for (j in i:dimension) {
      matrix[i,j] <- vec[i + j] # forces symmetry
      if ( matrix[i,j] < upper) {
        matrix[i,j] <- 0
      }
    }
  }
  diag(matrix) <- rep(1, dimension)
  # Now we ensure diagonal dominance
  for (k in 1:dimension) {
    matrix[k,] <- matrix[k,]/sum(abs(matrix[k,])) 
  }
  return(matrix)
}

set.seed(1)
n <- 100
precision <- matrix(c(4, 5, 0.01, # off-diagonal: s_xy <= s_x*s_y
                      5, 8, 0,
                      0.01, 0, 6), 3, 3)
mu0 <- c(2, 4, 8)
mat <- MASS::mvrnorm(n = n, mu = mu0, 
                     Sigma = solve(precision),
                     tol = 1e-8, empirical = TRUE)
lapply(c(1:3), function(i) eval(parse(text = paste0("x", i, " <<- mat[, ", i, "]"))))
y <- 100 - 4*x1 + 3*x2 - 2*x3 + rnorm(n, 0, 5)

df <- data.frame(x1 = x1, x2 = x2, x3 = x3, y = y, stringsAsFactors = FALSE)
plot(df)
par(mfrow = c(1, 2))
boxplot(df[, c(1:3)], names = c("x1", "x2", "x3"))
boxplot(df[, 4], xlab = "y")
par(mfrow = c(1, 1))
corrplot::corrplot(cor(df), type = "upper")

fit <- lm(formula = y ~ x1 + x2 + x3, data = df)
print(summary(fit))
# plenty of space for improvement :-)

In some way related Question: (1)

Answer 1

Here is a response inspired by Gaussian graphical models. Under gaussianity, multiple linear regression is intimately connected to the notion of partial correlation through the so-called precision matrix (inverse covariance matrix).

Let $\Omega$ denote the precision matrix of a set $\{X_1,...,X_r\}$ of variables. On the one hand, the coefficient $\beta_{j,k}$ of the regression of $X_j$ on $X_k$ is given by $-\frac{\Omega_{jk}}{\Omega_{jj}}$. (The diagonal elements $\Omega_{jj}$ are reciprocals of the conditional variances given the remaining variables.)

On the other hand, the $(r \times r)$ matrix $C$ whose $(j,k)$ entry equals $\frac{\Omega_{jk}}{\sqrt{\Omega_{jj}\Omega_{kk}}}$ is the negative partial correlation matrix (off the diagonal).

Now, if you can generate sparse positive definite matrices possibly with specific partial correlation, you simply use the mvrnorm function from the MASS package as follows:

mvrnorm(n=50,mu=rep(0,nrow(precision)),Sigma=solve(precision),tol=1e-8,empirical=TRUE),

where $precision$ denotes the precision matrix you generated.

Here is some code for generating a sparse definite positive matrix:

# input:
# dimension: number of rows of the precision matrix
# seed: random seed
# upper (in (0,1)) specifies the range of allowed non-zero entries 
# in the starting point for the construction of the precision matrix
generate.sparse.precision.matrix <- function(dimension,upper, seed){
  matrix <- matrix(rep(0,dimension*dimension), ncol=dimension)
  set.seed(seed)
  vec <- runif(n=dimension^2, min = 0, max = 1)

  for (i in 1:dimension){
    for (j in i:dimension){
      matrix[i,j] <- vec[i+j] # forces symmetry
      if( matrix[i,j] < upper){
        matrix[i,j] <- 0
      }
    }
  }
  diag(matrix) <- rep(1,dimension)
  # Now we ensure diagonal dominance
  for (k in 1:dimension){
    matrix[k,] <- matrix[k,]/sum(abs(matrix[k,])) 
  }
  return(matrix)
}

You can easily specify the partial correlations you want in the row(column) associated with the response variable in compliance with diagonal dominance.

EDIT:

Denote the standard deviation of the errors by $\sigma$. You have control over the variances of the predictors and the errors. Also, you know that the denominator in Student's statistic for coefficient $n°i$ equals $\sigma$ multiplied by the square root of the $i$th entry on the diagonal of the precision matrix of the predictors divided by the square root of the sample size.

You may for example play with the variances of the predictors given desired partial correlations and true underlying coefficients. This way, you will freely allow Student's statistic to vary for any given coefficient and it will naturally affect the p-values for the coefficients of the predictors.

As regards the p-value for the intercept, just shift the response variable by a parameter $t$ across a wide enough range of values. Remember that the intercept is just the mean of the response given that all predictors are equal to zero. Therefore, if you arrange for a model to be intercept-free and just re-iterate the estimations with different shifts of the response variable, then the larger $t$ (in absolute value) the higher the p-value.

I hope this helps.

EDIT (December): I just went over my answer after a while. Once you have generated data with your chosen partial correlations, you may simply replace each explanatory variable $X_i$ by $X_i+c_iX_i^2$ (or other exponents) with varying values for $c_i$. And then, you take any linear combination of your choosing as your response variable. The magnitudes of the $c_i's$ will affect the p-values. The larger in absolute value the more insignificant your new predictors $X_i+c_iX_i^2$ will be deemed.