2016-03-21 30 views
1

Ich habe eine Funktion geschrieben, um die Koeffizienten eines (linearen) Modells zu nehmen und sie auf die ursprünglichen Variablen anzuwenden, um einen Datenrahmen von Termen zu erhalten, die, wenn sie addiert werden, dem Ergebnis von vorhergesagt() entsprechen. Diese Fähigkeit scheint mir nützlich zu sein, um den Einfluss jeder Variablen (oder eines komplexeren Interaktionsterms usw.) auf das Modell besser zu verstehen.R: Modellkoeffizienten über Variablen im Modell anwenden, ein besserer Weg?

Gibt es einen besseren Weg? Ich fühle mich wie ein Hack. Ich habe mir die str() -Modelle angesehen und sehe bis jetzt keine einfachere Lösung. Der schwierige Teil besteht darin, Interaktionsterme zu erfassen und anzuwenden.

library(plyr) 
nospredict = function(model, data = model$model, sorted = TRUE) { # model is model (from lm, glm...), data is data.frame to be applied to                  
    c = coef(model) # model must support coef()                                        
    my.names = names(c) = gsub(':', '*', names(c)) # ':' equals multiplication in formulas, coefs                           
    data = data[ , colnames(data) %in% my.names] # don't do the attach() below with a zillion variables...                          
    final.out = adply(data, 1, function(y) { # did I mention I like plyr?                                  
     attach(as.list(y), warn.conflicts = FALSE) # so you can do eval algebra blackRmagic                             
     out = ldply(my.names, function (x) { # did I mention...                                    
      Intercept = 1 # (Intercept) from model is a constant, multiply it by 1                               
      eval( parse( text = paste(c[x], "*", x) ) )  }) # blackRmagic                               
     out = as.data.frame(t(out)) ; colnames(out) = my.names ; out 
    }) 
    rownames(final.out) = rownames(data) 
    final.out$Predict = predict(model, data) ## add predict() as column                                  
    if (sorted) { 
     final.out[order(final.out$Predict), ] ## return df sorted by predict()                                
    } 
    final.out 
} 
set.seed(10538) 
df = data.frame(a = 1:10, b = rnorm(10), c = 1:10 + rnorm(10)) 
lmf = lm(c ~ a * b, data = df) 

> df 
a   b   c 
1 1 -0.41184664 1.3739709 
2 2 1.06464586 0.8975101 
3 3 -0.07522363 3.4910425 
4 4 1.21643049 2.8856876 
5 5 0.34061917 4.3851439 
6 6 -1.00020786 6.1836535 
7 7 -0.36954963 6.4734150 
8 8 -1.47754640 6.8150569 
9 9 -0.19312147 9.6432687 
10 10 2.32220098 9.4276813 


> nospredict(lmf) 
(Intercept)   a   b   a*b Predict 
1 0.09801818 0.9282185 0.48332671 -0.05438652 1.4551769 
2 0.09801818 1.8564370 -1.24942570 0.28118420 0.9862137 
3 0.09801818 2.7846555 0.08827944 -0.02980103 2.9411521 
4 0.09801818 3.7128740 -1.42755405 0.64254425 3.0258824 
5 0.09801818 4.6410925 -0.39973700 0.22490279 4.5642765 
6 0.09801818 5.5693110 1.17380385 -0.79249635 6.0486367 
7 0.09801818 6.4975295 0.43368863 -0.34160685 6.6876294 
8 0.09801818 7.4257480 1.73398922 -1.56094237 7.6968130 
9 0.09801818 8.3539665 0.22663962 -0.22952439 8.4490999 
10 0.09801818 9.2821850 -2.72524198 3.06658890 9.7215501 
+0

letzte sortierung von vorhergesagt() werte hier nicht funktioniert, trivial fix, nicht wichtig, um die Frage. –

+0

Natürlich ist das reproduzierbare Beispiel oben nur das. Ich suche nach einer Lösung, die in den meisten oder allen Modellformeln verallgemeinerbar ist (ohne zu wissen, dass das Modell "a" und "b" enthält). –

Antwort

1
junk <- data.frame(x1 = rnorm(100), x2 = rnorm(100)) 

junk$YY <- 2 * junk$x1 + 4 * junk$x2 + 6 * junk$x1 * junk$x2 + 7 + rnorm(100) 

out <- lm(YY ~ x1 * x2, data = junk, x = TRUE) # The x = TRUE part is key! 

head(out$x) 
    (Intercept)   x1   x2  x1:x2 
1   1 -0.34885894 -0.8127228 0.283525629 
2   1 -0.04482498 -0.1601600 0.007179167 
3   1 -1.11721391 0.3266213 -0.364905892 
4   1 -0.08530188 0.3482372 -0.029705287 
5   1 0.19138684 -0.1659683 -0.031764149 
6   1 -1.89493717 1.0261454 -1.944481020 

coef(out) 
(Intercept)   x1   x2  x1:x2 
    7.053434 1.804441 4.130249 5.970553 

nomThings <- t(t(out$x[, names(coef(out))]) * coef(out)) 

Ich bin mir nicht ganz sicher, ob diese korrekt funktionieren, wenn Sie einige Faktoren als unabhängige Variablen haben, oder dass es in allen Situationen korrekt funktionieren. Aber ich vermute es.

Sie könnten natürlich coef (out) als temporäres Objekt speichern und so weiter, um den Code etwas lesbarer und etwas effizienter zu machen.

Angesichts, dass Sie dies leicht tun können, würde ich darüber nachdenken, ob Sie sollten.

+0

Danke. x = TRUE im Aufruf lm() ist die magische Option, die eine vollständige Modellmatrix bereitstellt und das Leben viel einfacher und wahrscheinlich zuverlässiger macht. Unabhängig davon, wie wichtig es ist, dies zu tun. –