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 <[email protected]> wrote:
> On Wed, Dec 22, 2010 at 12:39 PM, Advait Godbole
> <[email protected]> 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]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
