[R-sig-Geo] Breakpoint analysis with two variables

[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
     }
}


> z <- calc(z, fun)

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):

> d <-  z[10]
> fun(d)

Error in segmented.default(coef, ~s, psi = min(x[1:11], na.rm = TRUE)) :
  No default method for segmented

# or for a few cellls
> apply(d[1:5], 1, fun)
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:
> 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.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo