Skip to content

help with predict and plotting confidence intervals

4 messages · Michael Denslow, David Winsemius

Michael Denslow
# 
Dear R help,

This seems to be a commonly asked question and I am able to run examples that have been proposed, but I can't seems to get this to work with my own data. Reproducible code is below. Thank you in advance for any help you can provide.

The main problem is that I can not get the confidence lines to plot correctly.
The secondary problem is that predict is not able to find my object when 
I include a model object.

## THE DATA
wt.data <- data.frame(code = factor(LETTERS[1:24]),
	area = c(60865,480,656792,92298,1200,1490,8202,4000,220,245,4000,390,325,
	16,162911,20235,68800,3389,7,696,4050,1498,1214,99460),
	species = c(673,650,1353,1026,549,536,782,734,516,580,673,560,641,443,1105,
	871,789,575,216,407,942,655,582,1018))

# TRANSFORM AND ADD TO DATAFRAME
wt.data$logA <- log10(wt.data$area)
wt.data$logS <- log10(wt.data$species)

wt.mod <- lm(logS~logA, data = wt.data)

# PLOT THE DATA
with(wt.data,plot(logA,logS, ylim = c(2.3,3.2),xlim = c(0,6)))
abline(wt.mod, lwd = 2)


# create a prediction dataframe the same length as data
pred.frame <- data.frame(a = seq(0,6, length.out = 24))

# error ' object "logA" not found'
# I am not sure why object is not found, I assume this has to do with 
# the way I added the transformed variables to the dataframe
pp <- predict(wt.mod, int = "p", newdata=pred.frame)

# runs ok?
pp <- predict(lm(wt.data$logS~wt.data$logA), int = "p", newdata=pred.frame)

# lines are jagged??
# I am not sure how to get the lines to draw correctly here
matlines(pred.frame$a,pp, lty=c(1,2,2),col="black")
R version 2.8.1 (2008-12-22) 
i386-pc-mingw32 

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United States.1252;LC_MONETARY=English_United States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] tools_2.8.1

Michael Denslow

Graduate Student
I.W. Carpenter Jr. Herbarium [BOON]
Department of Biology
Appalachian State University
Boone, North Carolina U.S.A.

-- AND --

Communications Manager
Southeast Regional Network of Expertise and Collections
sernec.org
David Winsemius
#
On Mar 12, 2009, at 11:14 AM, Michael Denslow wrote:

            

      
      
      
    

        
It is commonly asked ... and commonly answered.

            

      
      
      
    

        
I suspect you omitted the actual call that produced this error. I  
suspect it was something along the lines of:
 > predwt <- predict(wt.mod, newdata=pred.frame)
Error in eval(expr, envir, enclos) : object "logA" not found

            

      
      
      
    

        
Because you didn't give the arguments the same name(s) as were used in  
the model formula.

            

      
      
      
    

        
Lines? What lines? If predict does not find a newdata object that  
satisfies its requirements, it uses the original data.

            

      
      
      
    

        
The x values are your sequence whereas the y values are in the  
sequence from the original data. They are not correctly associated  
with each other.

Try:
pp <- predict(lm(wt.data$logS~wt.data$logA), int = "p", newdata=  
data.frame(logA=seq(0,6, length.out = 24))  )
  plot(pp)

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      David Winsemius
    

    

      
      #
    

  


  

      
On Mar 12, 2009, at 11:45 AM, David Winsemius wrote:

        

            

      
      
      
    

        
At this point I should not have accepted your starting point. A better  
starting point would be to use the wt.mod model:

pp <- predict(wt.mod, int = "p", newdata= list(logA=seq(0,6,  
length.out = 24)) )

# Followed by:

plot( seq(0,6, length.out = 24),  pp[ ,"fit"] )
  lines(seq(0,6, length.out = 24),  pp[ ,"lwr"], lty=2)
  lines(seq(0,6, length.out = 24),  pp[ ,"upr"], lty=2)

            

      
      
      
    

        

      
      
    

        
David Winsemius, MD
Heritage Laboratories
West Hartford, CT

  
  


  


      

    

    

      
      

        


  

        

      Michael Denslow
    

    

      
      #
    

  


  

      
Thank you for you help Dr. Winsemius.

The problem seems to stem from the fact that I have used the incorrect name in the prediction dataframe.
The following code seems to work correctly. 
Thank you again,
Michael

wt.data <- data.frame(code = factor(LETTERS[1:24]),
	area = c(60865,480,656792,92298,1200,1490,8202,4000,220,245,4000,390,325,
	16,162911,20235,68800,3389,7,696,4050,1498,1214,99460),
	species = c(673,650,1353,1026,549,536,782,734,516,580,673,560,641,443,1105,
	871,789,575,216,407,942,655,582,1018))

wt.data$logA <- log10(wt.data$area)
wt.data$logS <- log10(wt.data$species)
wt.mod <- lm(logS~logA, data = wt.data)

with(wt.data,plot(logA,logS, ylim = c(2.0,3.5),xlim = c(0,6)))
pred.frame <- data.frame(logA = seq(0,6, length.out = 24))
pp <- predict(wt.mod, int = "p", newdata=pred.frame)
matlines(pred.frame$logA,pp, lty=c(1,2,2),col="red")




Michael Denslow

Graduate Student
I.W. Carpenter Jr. Herbarium [BOON]
Department of Biology
Appalachian State University
Boone, North Carolina U.S.A.

-- AND --

Communications Manager
Southeast Regional Network of Expertise and Collections
sernec.org