Как получить коэффициенты регрессии и соответствия модели, используя матрицу корреляции или ковариации вместо фрейма данных с использованием R?

Я хочу иметь возможность регрессировать коэффициенты из множественной линейной регрессии, предоставляя матрицу корреляции или ковариации вместо data.frame. Я понимаю, что вы теряете некоторую информацию, имеющую отношение к определению перехвата и т. Д., Но даже если матрица корреляции должна быть достаточной для получения стандартизированных коэффициентов и объясненных оценок дисперсии.

Так, например, если у вас были следующие данные

# get some data
library(MASS)
data("Cars93")
x <- Cars93[,c("EngineSize", "Horsepower", "RPM")]

Вы можете запустить регрессию следующим образом:

lm(EngineSize ~ Horsepower + RPM, x)

но что если вместо данных у вас есть корреляционная матрица или ковариационная матрица:

corx <- cor(x)
covx <- cov(x)
  • Какая функция в R позволяет вам запускать регрессию на основе матрицы корреляции или ковариации? В идеале это должно быть похоже на lm так что вы можете легко получить такие вещи, как г-квадрат, скорректированный г-квадрат, прогнозируемые значения и так далее. Предположительно, для некоторых из этих вещей вам также потребуется указать размер выборки и, возможно, вектор средних значений. Но это также было бы хорошо.

То есть что-то вроде:

lm(EngineSize ~ Horsepower + RPM, cov = covx) # obviously this doesn't work

Обратите внимание, что этот ответ на Stats.SE дает теоретическое объяснение того, почему это возможно, и предоставляет пример некоторого пользовательского кода R для вычисления коэффициентов?

4 ответа

Решение

Используя lavaan, вы можете сделать следующее:

library(MASS)
data("Cars93")
x <- Cars93[,c("EngineSize", "Horsepower", "RPM")]

lav.input<- cov(x)
lav.mean <- colMeans(x)

library(lavaan)
m1 <- 'EngineSize ~ Horsepower+RPM'
fit <- sem(m1, sample.cov = lav.input,sample.nobs = nrow(x), meanstructure = TRUE, sample.mean = lav.mean)
summary(fit, standardize=TRUE)

Результаты:

Regressions:
                   Estimate    Std.Err  Z-value  P(>|z|)   Std.lv    Std.all
  EngineSize ~                                                              
    Horsepower          0.015    0.001   19.889    0.000      0.015    0.753
    RPM                -0.001    0.000  -15.197    0.000     -0.001   -0.576

Intercepts:
                  Estimate    Std.Err  Z-value  P(>|z|)   Std.lv    Std.all
   EngineSize          5.805    0.362   16.022    0.000      5.805    5.627

Variances:
                  Estimate    Std.Err  Z-value  P(>|z|)   Std.lv    Std.all
    EngineSize          0.142    0.021    6.819    0.000      0.142    0.133

Я думаю, что Lavaan звучит как хороший вариант, отмечу, что @Philip указал мне правильное направление. Я только что упомянул здесь, как извлечь несколько дополнительных функций модели, используя lavaan (в частности, r-квадрат и скорректированный r-квадрат), которые вы можете захотеть.

Для получения последней версии см.: https://gist.github.com/jeromyanglim/9f766e030966eaa1241f10bd7d6e2812:

# get data
library(MASS)
data("Cars93")
x <- Cars93[,c("EngineSize", "Horsepower", "RPM")]

# define sample statistics 
covx <- cov(x)
n <- nrow(x)
means <- sapply(x, mean) # this is optional


fit <- lavaan::sem("EngineSize ~ Horsepower + RPM", sample.cov = covx,
                   sample.mean = means,
                    sample.nobs = n)

coef(fit) # unstandardised coefficients
standardizedSolution(fit) # Standardised coefficients
inspect(fit, 'r2') # r-squared

# adjusted r-squared
adjr2 <- function(rsquared, n, p) 1 - (1-rsquared)  * ((n-1)/(n-p-1))
# update p below with number of predictor variables
adjr2(inspect(fit, 'r2'), n = inspect(fit, "nobs"), p = 2) 

Пользовательская функция

А здесь есть небольшая функция, которая обеспечивает подгонку из лавы, а также несколько важных функций (т. Е. В основном упаковывает большинство из вышеперечисленных). Это предполагает случай, когда у вас нет средств.

covlm <- function(dv, ivs, n, cov) {
    # Assumes lavaan package
    # library(lavaan)
    # dv: charcter vector of length 1 with name of outcome variable
    # ivs: character vector of names of predictors
    # n: numeric vector of length 1: sample size
    # cov: covariance matrix where row and column names 
    #       correspond to dv and ivs
    # Return
    #      list with lavaan model fit
    #      and various other features of the model

    results <- list()
    eq <- paste(dv, "~", paste(ivs, collapse = " + "))
    results$fit <- lavaan::sem(eq, sample.cov = cov,
                       sample.nobs = n)

    # coefficients
    ufit <- parameterestimates(results$fit) 
    ufit <- ufit[ufit$op == "~", ]
    results$coef <- ufit$est
    names(results$coef) <- ufit$rhs

    sfit <- standardizedsolution(results$fit) 
    sfit <- sfit[sfit$op == "~", ]
    results$standardizedcoef <- sfit$est.std
    names(results$standardizedcoef) <- sfit$rhs

    # use unclass to not limit r2 to 3 decimals
     results$r.squared <- unclass(inspect(results$fit, 'r2')) # r-squared

    # adjusted r-squared
      adjr2 <- function(rsquared, n, p) 1 - (1-rsquared)  * ((n-1)/(n-p-1))
    results$adj.r.squared <- adjr2(unclass(inspect(results$fit, 'r2')), 
                                n = n, p = length(ivs)) 
    results

}

Например:

x <- Cars93[,c("EngineSize", "Horsepower", "RPM")]
covlm(dv = "EngineSize", ivs = c("Horsepower", "RPM"),
      n = nrow(x), cov = cov(x))

Это все производит:

$fit
lavaan (0.5-20) converged normally after  27 iterations

  Number of observations                            93

  Estimator                                         ML
  Minimum Function Test Statistic                0.000
  Degrees of freedom                                 0
  Minimum Function Value               0.0000000000000

$coef
 Horsepower         RPM 
 0.01491908 -0.00100051 

$standardizedcoef
Horsepower        RPM 
 0.7532350 -0.5755326 

$r.squared
EngineSize 
     0.867 

$adj.r.squared
EngineSize 
     0.864 

Помните, что:

$ Бета =(X'X)^-1. X'Y$

Пытаться:

(bs<-solve(covx[-1,-1],covx[-1,1]))

 Horsepower         RPM 
 0.01491908 -0.00100051 

Для перехвата вам понадобятся средние значения переменных. Например:

  ms=colMeans(x)
  (b0=ms[1]-bs%*%ms[-1])

         [,1]
[1,] 5.805301

Другим типичным решением является создание набора данных, который имеет ту же матрицу дисперсии-ковариации, что и исходные данные. Вы можете сделать это с mvrnorm() в MASS пакет. С помощью lm() для этого нового набора данных будут получены оценки параметров и стандартные ошибки, идентичные тем, которые были бы оценены из исходного набора данных (за исключением перехвата, который недоступен, если у вас нет средств для каждой переменной). Вот пример того, как это будет выглядеть:

#Assuming the variance covariance matrix is called VC
n <- 100 #sample size
nvar <- ncol(VC)
fake.data <- mvrnorm(n, mu = rep(0, nvar), sigma = VC, empirical = TRUE)
lm(Y~., data = fake.data)
Другие вопросы по тегам