Wie erhalten Sie Regressionskoeffizienten und Modellanpassungen unter Verwendung von Korrelations- oder Kovarianzmatrix anstelle von Datenrahmen unter Verwendung von R?

Question

Ich möchte in der Lage sein, Regressionskoeffizienten aus multipler linearer Regression durch Bereitstellen einer Korrelations- oder Kovarianzmatrix anstelle eines data.frame. Mir ist klar, dass Sie einige Informationen verlieren, die für die Bestimmung des Achsenabschnitts relevant sind, aber es sollte sogar die Korrelationsmatrix ausreichen, um standardisierte Koeffizienten und Varianzschätzwerte zu erhalten.

So zum Beispiel, wenn Sie die folgenden Daten haben

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

Sie eine Regression laufen könnten wie folgt aussehen:

lm(EngineSize ~ Horsepower + RPM, x) 

aber was, wenn statt Daten, die Sie die Korrelationsmatrix hatten oder die Kovarianzmatrix:

corx <- cor(x) 
covx <- cov(x) 

• Mit welcher Funktion in R können Sie eine Regression basierend auf der Korrelations- oder Kovarianzmatrix durchführen? Im Idealfall sollte es ähnlich sein wie lm, so dass Sie leicht Dinge wie r-Quadrat, angepasste r-Quadrat, vorhergesagte Werte und so weiter erhalten können. Vermutlich müssten Sie für einige dieser Dinge auch die Stichprobengröße und möglicherweise einen Vektor von Mitteln angeben. Aber das wäre auch in Ordnung.

Das heißt, so etwas wie:

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

Beachten Sie, dass diese Antwort auf Stats.SE eine theoretische Erklärung sieht, warum es möglich ist, und liefert ein Beispiel für einige Code benutzerdefinierte R für Koeffizienten Berechnung?

Answer 1

Mit lavaan Sie Folgendes tun könnte:

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) 

Ergebnisse sind:

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 

Answer 2

Denken Sie daran, dass:

$ beta = (X'X)^- 1. X'Y

$

Versuchen:

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

Horsepower   RPM 
0.01491908 -0.00100051 

Für die Intercept Sie werden Mittelwerte der Variablen benötigen. Zum Beispiel:

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

     [,1] 
[1,] 5.805301 

Answer 3

Ich denke lavaan wie eine gute Option klingt, stelle ich fest, dass @Philip wies mich in die richtige Richtung. Ich erwähne hier nur, wie man ein paar zusätzliche Modell-Features mit Hilfe von lavaan (insbesondere r-squared und adjusted r-squared) extrahieren kann, die Sie vielleicht wollen.

Die aktuelle Version finden Sie unter: 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) 

Benutzerdefinierte Funktion

Und hier ist ein bisschen eine Funktion, die den Sitz von lavaan zusammen mit einigen Funktionen von Bedeutung (dh liefert, im Grunde Verpackung der meisten der oben genannten).Es nimmt einen Fall an, in dem Sie nicht die Mittel haben.

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 

} 

Zum Beispiel:

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

Dies alles ergibt:

$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 

Answer 4

Eine andere Art von funky Lösung ist, einen Datensatz zu erzeugen, der den gleichen Varianz-Kovarianzmatrix wie das Original aufweist Daten. Sie können dies mit mvrnorm() im MASS Paket tun. Die Verwendung von lm() für diesen neuen Datensatz führt zu Parameterschätzungen und Standardfehlern, die mit denen identisch sind, die aus der ursprünglichen Datenmenge geschätzt worden wären (außer für den Achsenabschnitt, auf den nur zugegriffen werden kann, wenn die Mittelwerte jeder Variablen vorhanden sind). Hier ist ein Beispiel dafür, wie das aussehen würde:

#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)