2013-10-24 5 views
10

Angenommen, Sie modellieren binomische Daten, wobei jede Antwort eine Anzahl von Erfolgen (y) aus einer Anzahl von Versuchen (N) mit einigen erklärenden Variablen (a und b) ist. Es gibt ein paar Funktionen, die diese Art der Sache zu tun, und sie scheinen alle unterschiedliche Methoden zu verwenden, um anzugeben, y und N.Welche alternativen Möglichkeiten gibt es, binomiale Erfolge/Versuche in einer Formel zu spezifizieren?

In glm, tun Sie glm(cbind(y,N-y)~a+b, data = d) (Matrix des Erfolgs/Fehler auf LHS)

In InlA , tun Sie inla(y~a+b, Ntrials=d$N, data=d)

In glmmBUGS (Anzahl der Versuche separat angeben), tun Sie glmmBUGS(y+N~a+b,data=d) (angeben Erfolg + Studien als Begriffe auf LHS)

Wenn neue Methoden der Programmierung, ich habe immer gedacht, es am besten zu folgen, was Glm macht das, da würden die Leute normalerweise zuerst auf bin treffen omiale Antwortdaten. Ich kann mich jedoch nie daran erinnern, ob es cbind(y,N-y) oder cbind(y,N) ist - und ich habe normalerweise Erfolg/Anzahl der Versuche in meinen Daten anstatt Erfolg/Anzahl der Fehlschläge - YMMV.

Andere Ansätze sind natürlich möglich. Zum Beispiel eine Funktion auf der RHS mit markieren, ob die Variable Anzahl von Versuchen oder die Anzahl der fehlschlägt:

myblm(y ~ a + b + Ntrials(N), data=d) 
myblm(y ~ a + b + Nfails(M), data=d) # if your dataset has succ/fail variables 

oder definieren einen Operator nur einen cbind zu tun, so können Sie tun:

myblm(y %of% N ~ a + b, data=d) 

, wodurch der LHS eine gewisse Bedeutung verliehen wird.

Hat jemand bessere Ideen? Was ist der richtige Weg, dies zu tun?

+1

Ich mag% von%, weil es die Reihenfolge offensichtlich macht. Aber es ist weniger Standard .... –

+2

'glm' erlaubt auch die Proportion auf der LHS mit einem zusätzlichen' Gewicht' Argument geben die Nenner gegeben –

Antwort

-2

Von Hilfeseite auf GLM des r: „... oder als zweispaltige Matrix mit den Spalten, die Zahl der Erfolge und Ausfälle geben

So hat es cbind werden (Y , NY)

+2

Was fügt dies zur Frage des OP hinzu? Er hat das schon gesagt, denke ich. –

+0

Ich denke, es hat mit dem folgenden Kommentar des OP zu tun: _Ich kann mich nie erinnern, ob sein 'cbind (y, N-y)' oder 'cbind (y, N)' ..._. Es ist jedoch keine Antwort auf die Frage des OP. –

0

Ich mag diese Methode aus der GLM-Dokumentation:

Für binomische und quasibinomial Familien der Antwort auch als af angegeben werden kann Schauspieler (wenn die erste Ebene bezeichnet Versagen und alle anderen Erfolg)

Diese comports gut mit der Art und Weise Erfolgen und Misserfolgen oft in meiner Erfahrung entstehen. Einer ist ein Catch-All (z. B. "nicht gewählt") und es gibt eine Vielzahl von Möglichkeiten, um den anderen zu erreichen (z. B. "für A gewählt", "für B gewählt"). Ich hoffe, es ist klar aus der Art, wie ich das formuliere, dass "Erfolg" und "Versagen", wie definiert durch glm, willkürlich definiert werden kann, so dass die erste Ebene einem "Fehler" entspricht und alle anderen Ebenen ein "Erfolg" sind.

0

Sie können auch y Fraktion sein in diesem Fall müssen Sie die weights liefern. Es ist nicht in der formula Argument aber eine fast gleiche Anzahl von Tastenanschlägen, als ob es in der formula war.Hier ist ein Beispiel

> set.seed(73574836) 
> x <- runif(10) 
> n <- sample.int(10, 2) 
> y <- sapply(mapply(rbinom, size = 1, n, (1 + exp(1 - x))^-1), function(x) 
+ sum(x == 1)) 
> df <- data.frame(y = y, frac = y/n, x = x, weights = n) 
> df 
    y frac  x weights 
1 2 1.000 0.9051  2 
2 5 0.625 0.3999  8 
3 1 0.500 0.4649  2 
4 4 0.500 0.5558  8 
5 0 0.000 0.8932  2 
6 3 0.375 0.1825  8 
7 1 0.500 0.1879  2 
8 4 0.500 0.5041  8 
9 0 0.000 0.5070  2 
10 3 0.375 0.3379  8 
> 
> # the following two fits are identical 
> summary(glm(cbind(y, weights - y) ~ x, binomial(), df)) 

Call: 
glm(formula = cbind(y, weights - y) ~ x, family = binomial(), 
    data = df) 

Deviance Residuals: 
    Min  1Q Median  3Q  Max 
-1.731 -0.374 0.114 0.204 1.596 

Coefficients: 
      Estimate Std. Error z value Pr(>|z|) 
(Intercept) -0.416  0.722 -0.58  0.56 
x    0.588  1.522 0.39  0.70 

(Dispersion parameter for binomial family taken to be 1) 

    Null deviance: 9.5135 on 9 degrees of freedom 
Residual deviance: 9.3639 on 8 degrees of freedom 
AIC: 28.93 

Number of Fisher Scoring iterations: 3 

> summary(glm(frac ~ x, binomial(), df, weights = weights)) 

Call: 
glm(formula = frac ~ x, family = binomial(), data = df, weights = weights) 

Deviance Residuals: 
    Min  1Q Median  3Q  Max 
-1.731 -0.374 0.114 0.204 1.596 

Coefficients: 
      Estimate Std. Error z value Pr(>|z|) 
(Intercept) -0.416  0.722 -0.58  0.56 
x    0.588  1.522 0.39  0.70 

(Dispersion parameter for binomial family taken to be 1) 

    Null deviance: 9.5135 on 9 degrees of freedom 
Residual deviance: 9.3639 on 8 degrees of freedom 
AIC: 28.93 

Number of Fisher Scoring iterations: 3 

Der Grund, warum die oben genannten Arbeiten kommt darauf an, was glm tatsächlich tut für binomische Ergebnisse. Er berechnet einen Bruch für jede Beobachtung und ein Gewicht, das mit der Beobachtung verbunden ist, unabhängig davon, wie Sie das Ergebnis angeben. Hier ist ein Ausschnitt aus ?glm, die einen Hauch von gibt, was bei der Schätzung wird

Wenn ein binomialglm Modell, indem sie eine zweispaltige Antwort angegeben wurde, die zurück Gewichte von prior.weights sind die Gesamtzahl der Fälle (berücksichtigt durch die mitgelieferten Fallgewichte) und die Komponente y des Ergebnisses ist der Anteil der Erfolge.

Alternativ können Sie einen Wrapper machen für glm.fit oder glmmodel.frame verwenden. Siehe das ... Argument in ?model.frame

... für model.frame Methoden, eine Mischung aus weiteren Argumenten wie Daten, na.action, subset auf die Standardmethode zu übergeben. Etwaige zusätzliche Argumente (wie offset und weights oder andere benannte Argumente), die die Standardmethode erreicht, werden verwendet, um weitere Spalten im Modellrahmen, mit parenthesised Namen wie "(offset)".

Kommentar

Ich sah danach Ben Bolker Kommentar. Das oben genannte ist was er hervorhebt.

Verwandte Themen