Как получить коэффициенты регрессии и соответствия модели, используя матрицу корреляции или ковариации вместо фрейма данных с использованием 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)