2013-10-09 19 views
5

Ho un problema che sto cercando di risolvere senza successo. Più di due giorni di ricerche e non ho avuto un solo indizio. Scusate se la risposta è là fuori e non l'ho trovata.Come posso utilizzare la funzione di previsione in R in una regressione logistica installata anni fa?

Supponiamo di avere una regressione di equazioni logistica (modello binario) da un vecchio modello che avete stimato alcuni anni fa. Quindi conosci i parametri βk (k = 1, 2, ..., p) perché sono stati stimati in passato. Ma non hai i dati che sono stati usati per adattarsi al modello.

La mia domanda è: posso introdurre questo vecchio modello logistico stimato in R come oggetto (corrispondente a un modello di regressione logistica)?

Vorrei utilizzare la funzione "previsione" per provare questa regressione logistica con un nuovo set di dati (dati presenti) e quindi verificare la validità di questo vecchio modello in piedi nel tempo. E per utilizzare questa funzione è necessario l'oggetto del modello di regressione logistica.

Grazie mille in anticipo.

+0

Questa domanda sembra essere fuori tema, perché è sulle statistiche. Forse dovrebbe essere migrato a [Cross Validated] (http://stats.stackexchange.com). – Thomas

+7

L'utente sta provando a convertire un'equazione che (s) ha in un oggetto. Questa è una domanda abbastanza di programmazione, penso che si adatta bene. –

+0

Probabilmente modificherei un modello esistente, ma questo è imbroglio. –

risposta

6

Per il mio commento, penso che si potrebbe iniziare semplicemente calcolando le previsioni direttamente dai coefficienti. Ecco un esempio che mette a confronto l'output di predict.glm alla probabilità previste calcolate direttamente sui dati:

# construct some data and model it 
# y ~ x1 + x2 
set.seed(1) 
x1 <- runif(100) 
x2 <- runif(100) 
y <- rbinom(100,1,(x1+x2)/2) 
data1 <- data.frame(x1=x1,x2=x2,y=y) 
x3 <- runif(100) 
x4 <- runif(100) 
y2 <- rbinom(100,1,(x3+x4)/2) 
data2 <- data.frame(x1=x3,x2=x4,y=y2) 
glm1 <- glm(y~x1+x2,data=data1,family=binomial) 

# extract coefs 
#summary(glm1) 
coef1 <- coef(glm1) 

# calculate predicted probabilities for current data 
tmp1 <- coef1[1] + (data1$x1*coef1[2]) + (data1$x2*coef1[3]) 
pr1 <- 1/(1+(1/exp(tmp1))) 
# these match those from `predict`: 
all.equal(pr1,predict(glm1,data1,type='response')) 

# now apply to new data: 
tmp2 <- coef1[1] + (data2$x1*coef1[2]) + (data2$x2*coef1[3]) 
pr2 <- 1/(1+(1/exp(tmp2))) 
pr2 

Questo ovviamente non è una soluzione generale, né di gestire correttamente l'incertezza, ma penso che sia un approccio migliore rispetto l'hacking predict .

+0

Buona installazione - sarei tentato di eseguire 'qqplot' del set di dati rispetto ai dati" simulati "da questo modello. –

5

È possibile creare un glm fit con solo un offset creato dai coefficienti di cui si dispone, quindi utilizzare la funzione di previsione normale con quella. Ad esempio, utilizzando i dati dell'iride (primo montaggio di un modello sui dati reali, quindi il montaggio di un nuovo modello utilizzando dati fittizi ed i coefficienti della prima forma):

fit1 <- glm(I(Species=='versicolor') ~ Petal.Length + Petal.Width, 
    data=iris, family=binomial) 
coef(fit1) 

dummydata <- data.frame(Petal.Length = rnorm(10), Petal.Width=rnorm(10), 
    Species = rep(c('versicolor','other'), each=5)) 

fit2 <- glm(I(Species=='versicolor') ~ 0 + 
    offset(-2.863708 + 1.563076*Petal.Length - 3.153165*Petal.Width), 
    data=dummydata, family=binomial) 

pred1 <- predict(fit1, newdata=iris) 
pred2 <- predict(fit2, newdata=iris) 
plot(pred1,pred2) 
abline(0,1, col='green')