Skip to content

Breakpoint analysis with two variables

4 messages · nate_m, Robert J. Hijmans

Original
nate_m
# 
Hi all,

I am a beginner, so bear with me please. My goal is to perform breakpoint
analysis (package 'segmented') on each cell of two raster stacks (11 files
thick each). Some cells throughout both bricks are NA, so I have to skip
these as well. My goal is to output breakpoints with confidence intervals
based on the two variables for each cell. Here is my code and the error, I
think I am doing something incorrectly with data types;

library(raster)
#each stack is 11 rasters thick
y1 <- stack(y2)
x1 <- stack(x2)

#testing on small region
y <- crop(y1, extent(-103,-102,37,38))
x <- crop(x1, extent(-103,-102,37,38))

fun1 <- function(x) {
  m <- NA
  try(m <- lm(x ~ y), silent=T)
  m
}

library(segmented)

fun2 <- function(x) {
  m <- NA
  try(m <- segmented(x, ~y, psi=150), silent=T)
  m
}

fun3 <- function(x) {
  m <- NA
  try(m <- confint(x), silent=T)
  m
}

e1 <- calc(y, fun1)
os <- calc(e1, fun2)
results <- calc(os, fun3)

This is the error I receive after running the first function with e1;

Error in model.frame.default(formula = x ~ y, drop.unused.levels = TRUE) : 
  invalid type (S4) for variable 'y'

If anybody has a thought as to what I am doing wrong, it would be of great
help! Additionally, I would like to set "psi" from the segmented package
with values from a raster- is this possible?

Thanks in advance,
Nate



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Breakpoint-analysis-with-two-variables-tp7585607.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Robert J. Hijmans
# 
Nate,
With questions like this, please provide a fully reproducible example
with data created by the script or from example data from a package
(see the examples in the help files and the list archives for how to
that).

The error occurs because fun1 returns either NA or an object of class
'lm' (instead of a number). See ?calc for how to use lm with calc.
Robert
On Wed, Jan 22, 2014 at 12:48 PM, nate_m <Nathaniel.L.Mikle-1 at ou.edu> wrote:
5 days later
nate_m
# 
Hi Robert,

Thanks for taking the time to provide suggestions. Below is what I've worked
out, and I can now retrieve variables from lm (in this case the two
coefficients), but can't figure out how to save an entire "lm object" for
each cell to feed into the "segmented" portion of my script. So basically,
I'm looking to somehow save an entire "lm object" in raster brick form to
then utilize in the next step (where library(segmented) begins). Any ideas?

library(raster)
x <- raster(nrow=10, ncol=10)
y <- raster(nrow=10, ncol=10)
r <- stack( sapply(1:11, function(i) setValues(x, rnorm(ncell(x), i, 3) )) )
s <- stack( sapply(1:11, function(i) setValues(y, rnorm(ncell(y), i, 3) )) )
z <- stack(r,s)
z[1] <- NA
fun <- function(x) { if (is.na(x)) { return(cbind(NA,NA)) } else lm(x[1:11]
~ x[12:22])$coefficients }

g2 <- calc(z, fun)
g2

library(segmented)
min_x <- calc(r, min)
fun2 <- function(x) { if (is.na(x)) { return(cbind(NA,NA,NA,NA)) } else
segmented(x, ~s, psi = min_x)$coefficients }

seg_g2 <- calc(g2, fun2)



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Breakpoint-analysis-with-two-variables-tp7585607p7585654.html
Sent from the R-sig-geo mailing list archive at Nabble.com.
Robert J. Hijmans
# 
Nate,

I do not think you need to save the lm object. Instead you can write a
function that does the whole process. It appears your intention is
something like this:

fun <- function(x) {
    if ( is.na(x[1])  ) {  # or a variation on that such as  all(is.na(x))
         return( c(NA, NA,NA, NA) )
     } else {
         coef <- lm(x[1:11] ~ x[12:22])$coefficients
         segmented(coef, ~s, psi = min(x[1:11], na.rm=TRUE))$coefficients
     }
}
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
  cannot use this function

# That does not work, but that is probably because your call to
segmented in incorrect. The function should normally work for a single
cell (and that is how you can test it):
Error in segmented.default(coef, ~s, psi = min(x[1:11], na.rm = TRUE)) :
  No default method for segmented

# or for a few cellls
Error in apply(d[1:5], 1, fun) : dim(X) must have a positive length
Robert
On Tue, Jan 28, 2014 at 12:20 PM, nate_m <Nathaniel.L.Mikle-1 at ou.edu> wrote: