Dear Hadley,

Starting this in a new thread since I had tried something to that effect
like so:

After:

datCoef = adply( as.array(s_crop) , 1:2 ,
              function(x ) {
                            resp = x[1:120]
                            expl = x[121:240]
                            mdl = rlm(resp ~ expl)
                            res = c(mdl$coef[1] , mdl$coef[2])
                            #plt1 = plot(newx , predY)
                            #nam = x[ ,1]+x[ ,2]
                            #pdf(plt1 , paste(paste("plot",nam,
sep="_"),".pdf",sep="") )
                            return(res)

                          }
                          )

predY = rowSums(datCoef[ , 3:4] * as.matrix( rep(1, times = 4*4) , new_X ))
##  newx is the "future" x  vector)

Here I was trying to get the new Y by performing the operation "a + bX"
where X is the new 4 x 4 map and a (the intercept) and b (slope) should both
be 4 x 4

thanks,

Advait

On Wed, Dec 22, 2010 at 8:36 PM, Hadley Wickham <hadley@rice.edu> wrote:

> On Wed, Dec 22, 2010 at 12:39 PM, Advait Godbole
> <advaitgodbole@gmail.com> wrote:
> > Dear Members,
> >
> > I was able to find a solution after some deliberation and posting it here
> > for the benefit of forum members.
> >
> > To recap, my data was in the form of 2 netCDF files, each of dimension
> 356 x
> > 507 x 120 (lat,lon,time). I wanted to perform regressions between the two
> > variables on each grid point *over time*. I used "brick" to read in the
> > netcdfs as raster layers and stack along with the "*plyr"* package to
> > perform these regressions on the grid.
> >
> > # create raster object from the netcdf
> > rcmraster <- brick("regcm.cumul.hdd.10k.96_05.fixed.latlonvars.nc",
> varname
> > = "hdd")
> >
> > emitraster <- brick("
> > flip.percap.vulcan.10k.dp.build.381.RES.EIAmy.96_05.this.latest.nc
> > ",varname="emitpercap")
> >
> > # crop subset of above rasters, test method on subset which is 4 x 4 x
> 120
> > lim <- extent(rcmraster,180,185,300,305)
> > rcmraster_crop <- crop(rcmraster,lim)
> > emitraster_crop <- crop(emitraster,lim)
> >
> > # do some plotting to check
> >  par(mfrow=c(1,2))
> >  plot(rcmraster_crop[[1]])
> >  plot(emitraster_crop[[1]])
> >
> > # test regressions for cropped subset
> > s_crop <- stack(emitraster_crop,rcmraster_crop)
> >
> > # use of plyr (provides a more elegant alternative to lapply)
> > # the data frame datCoef contains the slope and intercept for each point
> in
> > the grid (denoted by X1 and X2)
> > datCoef = adply( as.array(s_crop) , 1:2 ,
> >              function(x ) {
> >                            resp = x[1:120]
> >                            expl = x[121:240]
> >                            mdl = lm(resp ~ expl)
> >                            res = c(mdl$coef[1] , mdl$coef[2])
> >                            return(res)
> >
> >                          }
> >                          )
>
> I'd do this in two steps:
>
> models <- alply( as.array(s_crop) , 1:2, function(x) {
>   resp = x[1:120]
>  expl = x[121:240]
>   lm(resp ~ expl)
> }
>
> Then you can access coefficients with
>
>  ldply(models, coef)
>
> and predictions with
>
>  laply(models, predict)
>
> but I think you might need to manipulate the output of predict to get
> the matrix shape that you are expecting.
>
> Hadley
>
> --
> Assistant Professor / Dobelman Family Junior Chair
> Department of Statistics / Rice University
> http://had.co.nz/
>



-- 
advait godbole

	[[alternative HTML version deleted]]

