[R-sig-Geo] Raster regression cooeficients plot?

[Robert Hijmans]

Hi Eddie,

Consider this:

> set.seed(0)
> y = runif(365)
> time = 1:365
> m = lm(y~time)
> m$coefficients
(Intercept)          time
 0.5297172419 -0.0001886475

But you do:

> summary(m)$coefficients

                 Estimate   Std. Error   t value     Pr(>|t|)
(Intercept)  0.5297172419 0.0292052249 18.137756 8.474764e-53
time        -0.0001886475 0.0001383047 -1.363999 1.734129e-01

Hence the 8 layers
> as.vector(summary(m)$coefficients)
[1]  5.297172e-01 -1.886475e-04  2.920522e-02  1.383047e-04
1.813776e+01 -1.363999e+00  8.474764e-53  1.734129e-01

I think you want to use this function (like the one in your original message!)

fun = function(x) { if (is.na(x[1])) { c(NA,NA) } else { lm(x ~
time)$coefficients } }

Note the two NA values such that the function always returns the same
number of values.
Before you use it in calc, you can first test it to see if it returns
what you expect

> fun(runif(365))
  (Intercept)          time
 0.5306391079 -0.0001631673


> funold=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time);summary(m)$coefficients}}
> funold(runif(365))
                 Estimate  Std. Error    t value     Pr(>|t|)
(Intercept)  4.814148e-01 0.030698435 15.6820636 1.099704e-42
time        -3.910719e-05 0.000145376 -0.2690073 7.880769e-01
>

Best, Robert


On Thu, Jun 16, 2011 at 1:29 AM, eddie [via R-sig-geo]
<ml-node+6482176-1417003307-149542 at n2.nabble.com> wrote:
> Dear Robert,
>
> Thank you very much for your reply. Sorry for not giving you enough
> information in order for you to give comments on what are the 8 coefficients
> generated. The full code is as below. My rasterbrick contains 365 layers.
>
> library(raster)
> setwd("C:/mydata")
> getwd()
> # define raster layer
> r <- raster(ncol=69, nrow=43, xmn=96.5, xmx=165.5, ymn=-16.5, ymx=26.5)
> # build rasterstack and then rasterbrick
> s = stack(list.files(pattern='*.rst'))
> b = brick(s)
> # fit linear model
> time <- 1:nlayers(b)
> fun=function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time);
> summary(m)$coefficients}}
> x3 <- calc(b, fun)
> plote(x3)
> # end
>
> Thank you.
>
>
> ________________________________
> If you reply to this email, your message will be added to the discussion
> below:
> http://r-sig-geo.2731867.n2.nabble.com/Raster-regression-cooeficients-plot-tp6478377p6482176.html
> To unsubscribe from Raster regression cooeficients plot?, click here.


--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Raster-regression-cooeficients-plot-tp6478377p6482308.html
Sent from the R-sig-geo mailing list archive at Nabble.com.