2016-06-22 21 views
1

Ich versuche, ein Oberflächendiagramm mit überlagerten Punkten aus einem Binomial-GLM mit der Funktion scatter3D() zu erstellen.GLM-Vorhersage für Oberflächendiagramm in Scatter3D() in R

Um dies zu tun, verwende ich Predic(), um die Z-Oberfläche für verschiedene Werte von x und y vorherzusagen.

# Data: 

library(plot3D) 

structure(list(
x = c(0.572082281112671, -0.295024245977402, 0.295024245977402, 0.861117839813232, 0.572082281112671, -1.74020183086395, 0.861117839813232, 0.283046782016754, 0.861117839813232, 0.283046782016754, -0.295024245977402, 1.43918883800507, 1.43918883800507, -0.295024245977402, -0.00598874036222696, -0.873095273971558, -0.295024245977402, -0.00598874036222696, -0.00598874036222696, 0.861117839813232), 
y = c(-1.09869265556335, -1.18406093120575, -0.0542464517056942, -0.192688703536987, -0.0208134315907955, 0.194501429796219, -0.126082852482796, 0.861439049243927, 0.624606966972351, -0.227061957120895, -1.32208430767059, -0.553429543972015, 0.538678884506226, 1.53797924518585, 0.230196505784988, 0.2959825694561, 0.158534705638885, 1.33240795135498, 0.0964559689164162, 0.740677952766418), 
z = c(0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 
w = structure(c(2L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L), .Label = c("0", "1"), class = "factor")), .Names = c("x", "y", "z", "w"), row.names = c(NA, 20L), class = "data.frame") 

Modell usw.

fit <-glm(formula = z ~ x * y + w, family = binomial) 

# x is continuous 
# y is continuous 
# w is dichotomous (yes, no, i.e. 0,1) [but see solution below] 
# z is dichotomous, but kept as numeric for plotting 

grid.lines = 100 
x.pred <- seq(min(x), max(x), length.out = grid.lines) 
y.pred <- seq(min(y), max(y), length.out = grid.lines) 
xy <- expand.grid(x = x.pred, y = y.pred) 

z.pred <- matrix(exp(predict(fit, newdata = xy)), 
      nrow = grid.lines, ncol = grid.lines) 

# fitted points for droplines to surface 
fitpoints <- exp(predict(fit)) 

Allerdings erhalte ich diese Fehlermeldung:

Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels) : variable lengths differ (found for 'w') 

W ist eine dritte Variable, die wichtig ist im Modell zu halten, aber ich kann‘ t finde heraus, wie man es konstant hält, während man noch die anderen Variablen aufträgt. Ich verstehe, dass ich etwas zwicken muss, aber ich kann nicht genau herausfinden, was das ist.

Beachten Sie, dass ich die Werte potenzieren, so dass sie eine Skala sind, die zwischen 0 und 1 eine Wahrscheinlichkeit ergibt, wenn ich sie grafisch darstelle. Wenn das nicht stimmt, lass es mich wissen. [Das war falsch - wies unten in den Kommentaren out]

ich damit fertig:

scatter3D(x, y, z, pch = 21, type = "p",col=rgb(red=0, green=17, blue=255, maxColorValue = 255, alpha = 150), bg = "#FF0000", 
     ylab = "Z-AM-Testosterone", xlab = "Z-AR-CAGn", zlab = "Divorce", 
     theta = -70, phi = 20, ticktype = "detailed", 
     surf = list(x = x.pred, y = y.pred, z = z.pred, 
        fit = fitpoints)) 

Ich bin sicher, dass es einfach ist, aber wenn jemand erklären könnte, wie w von der Vorhersage zu entfernen oder sie halten konstant, damit ich weitermachen kann, wäre ich sehr verpflichtet. Bitte schlagen Sie keine andere 3D-Methode vor - scatter3D ist besser als visreg oder andere für meine Zwecke.

Vielen Dank im Voraus für Ihre Hilfe.

+0

Können Sie bitte Daten enthalten, die uns mit einem [reproduzierbaren Beispiel] liefert (http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)? –

+1

Was passiert, wenn Sie einfach 'xy $ w <- ' (wobei "typisch" der Mittelwert, Median oder ein anderer vernünftiger Grundlinienwert ist) vor dem Vorhersagen verwenden? –

+0

Entschuldigung, lange Zeit lurker neues Plakat. Beispieldaten zur Verfügung gestellt, hoffentlich wie erforderlich (mein tatsächlicher Datensatz N = 675). – CPR

Antwort

1

Danke für die einfache Lösung, @Ben Bolker.

nahm ich den Mittelwert der numerischen Äquivalent der ja/nein, 0-1 Variable und geplottet nur die Vorhersagen aus, dass:

xy <- expand.grid(x = x.pred, y = y.pred, w = mean(w)) 

Dies erlaubte mir eine anständige aussehende Grafik zu erzeugen, die Sinn macht gegeben die Daten, unten gezeigt.

Scatter3D für negative Binomialmodell mit dem Mittelwert einer dritten dichotomous Kovariable (w), nach w in numerische Umwandlung:

3D plot of interaction