Plotting R squared in a loop

Question

My original dataset "d" includes monthly ndvi and precipitation data for 13 years for 26 stations:

head(d)

     row.names timestamp  station   year month   ndvi  landcover    altitude precipitation
1         1         1  A            2000   jan 0.4138 Mixed forest     2143          16.0
2      1769         2  A            2000   feb 0.4396 Mixed forest     2143           4.0
3      2055         3  A            2000   mar 0.4393 Mixed forest     2143          25.5
4      2341         4  A            2000   apr 0.6029 Mixed forest     2143          72.6
5      2627         5  A            2000   may 0.4756 Mixed forest     2143         241.7
6      2913         6  A            2000   jun 0.4969 Mixed forest     2143         505.9

The code below generates 6 plots that show the correlation of jan-jun precip to august ndvi. Now, I would like to add the r squared value to each plot, however doing the way I did in the code does not seem to work and results in r2=NA for each plot. I tried to get the r2 the same way as the p value but that also resulted in NAs.

I did the same thing for one station and the r2 worked. Does anyone have any ideas as to why it doesn't work for the code below?

d <- read.csv("a.csv", header = TRUE, sep = ",")
d <- na.omit(d)

for(m in c("jan","feb","mar","apr","may","jun")) {
  ndvi<-d$ndvi[d$month=="aug"]
  precip<-d$precipitation[d$month==m]
  r2<-cor(ndvi,precip)^2
  cat("month =",m,"P=",cor.test(ndvi,precip)$p.value,"\n")
  plot(ndvi~precip,main=m,
     sub=sprintf("r2=%.2f",r2)) 
  abline(lm(ndvi~precip))
}

The error I'm getting is:

Error in cor(ndvi, precip) : incompatible dimensions

Thanks for any help!

Answer 1

The problem is the NAs in the data - you need to keep them in until you've paired up the observations, but themn tell cor() and cor.test() what to do with the NA data - slightly confusingly they both have different ways of specifying that the NAs should be removed.

Does this work for your dataset now?

d <- read.csv("a.csv", header = TRUE, sep = ",")
#d <- na.omit(d) #keep NAs for the moment

for(m in c("jan","feb","mar","apr","may","jun")) {
  ndvi<-d$ndvi[d$month=="aug"]
  precip<-d$precipitation[d$month==m]
  r2<-cor(ndvi,precip,use="complete.obs")^2
  cat("month =",m,"P=",cor.test(ndvi,precip,na.action=na.omit)$p.value,"\n")
  plot(ndvi~precip,main=m,
     sub=sprintf("r2=%.2f",r2)) 
  abline(lm(ndvi~precip))
}