Podczas symulacji danych wielowymiarowych dla regresji, jak ustawić R-kwadrat (przykładowy kod w zestawie)?

Question

Próbuję symulować zestaw trzech zmiennych, aby móc na nim uruchomić modele regresji liniowej. "X1" i "X2" byłyby ciągłymi zmiennymi niezależnymi (średnia = 0, sd = 1), a "Y" byłaby ciągłą zmienną zależną.

Zmienne jest model regresji spowoduje współczynników tak: Y = 5 + 3 (X 1) - 2 (X2)

ja jak symulować zestawu danych, tak, że uzyskany model regresji ma R -dana wartość 0.2. Jak mogę określić wartość "sd.value", aby model regresji miał ten R-kwadrat?

n <- 200 
set.seed(101) 
sd.value <- 1 

X1 <- rnorm(n, 0, 1) 
X2 <- rnorm(n, 0, 1) 
Y <- rnorm(n, (5 + 3*X1 - 2*X2), sd.value) 

simdata <- data.frame(X1, X2, Y) 

summary(lm(Y ~ X1 + X2, data=simdata)) 

Answer 1

Spójrz na ten kod, powinien on być na tyle blisko, aby to, co chcesz:

simulate <- function(n.obs=10^4, beta=c(5, 3, -2), R.sq=0.8) { 
    stopifnot(length(beta) == 3) 
    df <- data.frame(x1=rnorm(n.obs), x2=rnorm(n.obs)) # x1 and x2 are independent 
    var.epsilon <- (beta[2]^2 + beta[3]^2) * (1 - R.sq)/R.sq 
    stopifnot(var.epsilon > 0) 
    df$epsilon <- rnorm(n.obs, sd=sqrt(var.epsilon)) 
    df$y <- with(df, beta[1] + beta[2]*x1 + beta[3]*x2 + epsilon) 
    return(df) 
} 
get.R.sq <- function(desired) { 
    model <- lm(y ~ x1 + x2, data=simulate(R.sq=desired)) 
    return(summary(model)$r.squared) 
} 
df <- data.frame(desired.R.sq=seq(from=0.05, to=0.95, by=0.05)) 
df$actual.R.sq <- sapply(df$desired.R.sq, FUN=get.R.sq) 
plot(df) 
abline(a=0, b=1, col="red", lty=2) 

Zasadniczo Twoje pytanie sprowadza się do zastanawianie się wyrażenie var.epsilon. Ponieważ mamy y = b1 + b2 * x1 + b3 * x2 + epsilon, a Xs i epsilon są niezależne, mamy var [y] = b2^2 * var [x1] + b3^2 * var [x2] + var [eps], gdzie var [Xs] = 1 z założenia. Możesz wtedy rozwiązać dla var [eps] jako funkcję R-kwadrat.

Answer 2

więc preparat do R^2 oznacza grupę 1-var (resztkowego)/var (łącznie)

W tym przypadku, wariancja Y będzie 3^2+2^2+sd.value^2, ponieważ dodawania trzech niezależnych zmiennych losowych. I asymptotycznie, wariancja resztkowa będzie po prostu sd.value^2.

Więc można obliczyć rsquared wyraźnie z tej funkcji:

rsq<-function(x){1-x^2/(9+ 4+x^2)} 

Przy odrobinie algebry można obliczyć odwrotność tej funkcji:

rsqi<-function(x){sqrt(13)*sqrt((1-x)/x)} 

Więc ustawienie sd.value<-rsqi(rsquared) powinno dać to, czego chcieć.

Możemy to sprawdzić w następujący sposób:

simrsq<-function(x){ 
    Y <- rnorm(n, (5 + 3*X1 - 2*X2), rsqi(x)) 
    simdata <- data.frame(X1, X2, Y) 
    summary(lm(Y ~ X1 + X2, data=simdata))$r.squared 
} 

> meanrsq<-rep(0,9) 
> for(i in 1:50) 
+ meanrsq<-meanrsq+Vectorize(simrsq)((1:9)/10) 
> meanrsq/50 
[1] 0.1031827 0.2075984 0.3063701 0.3977051 0.5052408 0.6024988 0.6947790 
[8] 0.7999349 0.8977187 

Wygląda więc na to, aby być poprawne.

Answer 3

To w jaki sposób to zrobić (ślepego iteracyjny algorytm, zakładając, że nie wiedzą, bo kiedy jesteś czysto zainteresowany „jak symulować to”):

simulate.sd <- function(nsim=10, n=200, seed=101, tol=0.01) { 
    set.seed(seed) 
    sd.value <- 1 
    rsquare <- 1:nsim 
    results <- 1:nsim 
    for (i in 1:nsim) { 
    # tracking iteration: if we miss the value, abort at sd.value > 7. 
    iter <- 0 
    while (rsquare[i] > (0.20 + tol) | rsquare[i] < (0.2 - tol)) { 
     sd.value <- sd.value + 0.01 
     rsquare[i] <- simulate.sd.iter(sd.value, n) 
     iter <- iter + 1 
     if (iter > 3000) { break } 
    } 
    results[i] <- sd.value # store the current sd.value that is OK! 
    sd.value <- 1 
    } 
    cbind(results, rsquare) 
} 

simulate.sd.iter <- function(sd.value, n=200) { # helper function 
    # Takes the sd.value, creates data, and returns the r-squared 
    X1 <- rnorm(n, 0, 1) 
    X2 <- rnorm(n, 0, 1) 
    Y <- rnorm(n, (5 + 3*X1 - 2*X2), sd.value) 
    simdata <- data.frame(X1, X2, Y) 
    return(summary(lm(Y ~ X1 + X2, data=simdata))$r.squared) 
} 

simulate.sd() 

kilka rzeczy do uwaga:

• Pozwoliłem, aby X1 i X2 były różne, ponieważ wpływa to na szukaną sd.value.
• Tolerancja określa dokładność tego oszacowania. Wszystko dobrze z r-kwadratem ~ 0,19 lub ~ 0,21? Czy tolerancja wynosi 0,01.
• Należy pamiętać, że zbyt precyzyjna tolerancja może uniemożliwić znalezienie wyniku.
• Wartość 1 jest dość złą wartością początkową, dzięki czemu ten iteracyjny algorytm jest dość powolny.

Otrzymany wektor 10 wyników jest:

[1] 5.64 5.35 5.46 5.42 5.79 5.39 5.64 5.62 4.70 5.55,

które trwa około 13 sekund w moim maszyny.

Mój następny krok to zacząć od 4,5, dodać 0,001 do iteracji zamiast 0,01 i być może obniżyć tolerancję. Powodzenia!

porządku, niektóre statystyki podsumowujące dla nsim = 100, biorąc 150 sekund z etapów wzrostu 0.001 i tolerancji jeszcze na 0,01:

Min. 1st Qu. Median Mean 3rd Qu. Max. 
4.513 4.913 5.036 5.018 5.157 5.393 

Dlaczego jesteś zainteresowany to chociaż?

Answer 4

Oto kolejny kod, aby wygenerować wielokrotnej regresji liniowej z błędów śledzić rozkład normalny: OPS Niestety ten kod tylko produkuje wielokrotnej regresji

sim.regression<-function(n.obs=10,coefficients=runif(10,-5,5),s.deviation=.1){ 
 
    
 
    n.var=length(coefficients) 
 
    M=matrix(0,ncol=n.var,nrow=n.obs) 
 
    
 
    beta=as.matrix(coefficients) 
 
    
 
    for (i in 1:n.var){ 
 
    M[,i]=rnorm(n.obs,0,1) 
 
    } 
 
    
 
    y=M %*% beta + rnorm(n.obs,0,s.deviation) 
 
    
 
    return (list(x=M,y=y,coeff=coefficients)) 
 
    
 
}