[R] Using Predict and GLM

[Bill.Venables at csiro.au]

Your problem is that your newdata data frame only specifies what the
value for A is in the predictions.  The values for W1 and W2 are
unspecified.  To predict from fitted models, you need to specify the
values you wish to use for *all* the predictor variables, not just the
one (I presume) is different from the estimation set.

Assuming that is what you want to do, here is a version of your script
that you may care to look at fairly closely.  I have omitted all the
unnecessary bits, (I could find), put the variables in a data frame
instead of leaving them strewn around the global environment, and done a
couple of predictions just to show it works.

_________________________________
set.seed(644)
n0 <- 200 # number of observations

myData <- data.frame(W1 = rnorm(n0, mean = 2, sd = 2),
  W2 = rnorm(n0, mean = 3, sd = 8))
myData <- transform(myData, # add A
  A = rbinom(n0, 1, ifelse(W1 > 1.5, 0.4, 0.2)))
myData <- transform(myData, # add Y
  Y = rbinom(n0, 1, 1/(1+exp(-(10*A-5*(W1)^2+2*W2)))))

Q <- glm(Y ~ A + W1 + W2, family = binomial, data = myData)

pData <- transform(myData, A = 0) # keep W1 and W2 as they were

QA <- predict(Q, newdata = pData) # linear predictors
QR <- predict(Q, newdata = pData, type = "response") # probs

plot(QA, QR)
_________________________________

Bill Venables. 


Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                         +61 4 8819 4402
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of Sherri Rose
Sent: Tuesday, 29 January 2008 1:27 PM
To: r-help at r-project.org
Subject: [R] Using Predict and GLM

Dear R Help,

I read through the archives pretty extensively before sending this  
email, as it seemed there were several threads on using predict with  
GLM.  However, while my issue is similar to previous posts (cannot  
get it to predict using new data), none of the suggested fixes are  
working.

The important bits of my code:

set.seed(644)
n0=200 #number of observations
W1=rnorm(n0,mean=2,sd=2) #Use rnorm to generate W1
W2=rnorm(n0,mean=3,sd=8) #Use rnorm to generate W1
Aprob=matrix(.2, nrow=n0, ncol=1) #generating the probability of A
#generating probability of A dependant on W1
for(i in 1:n0){
if (W1[i]>1.5) {Aprob[i]=0.4}
}
A=matrix(rbinom(n0, 1, Aprob), nrow=n0, ncol=1) #generating the 0/1  
exposure
Yprob=1/(1+exp(-(10*A-5*(W1)^2+2*W2)))
Y=matrix(rbinom(n0, 1, Yprob), nrow=n0, ncol=1) #generating the 0/1  
exposure
zero=data.frame(rep(0, n0))


Q=glm(cbind(Y, 1-Y) ~ A + W1 + W2, family='binomial')
QA=predict(Q, newdata=as.data.frame(A))
Q0=predict(Q,newdata=(A=zero))

I've tried many variations of the last line (Q0) to get the predicted  
values when A=0 with no luck.  With this code,  I get errors that my  
A=zero is a list even though I made it into a data frame.  This is  
the version of the code (after my reading) that *should* work for  
predict once I can get it to accept that it is not a list.

With other variants of the line that will run but are not  
syntactically correct, my QA and Q0 are the same.

Any guidance would be appreciated!

Sherri




	[[alternative HTML version deleted]]

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.