get.residuals.RdComputes (scaled) residuals of a set of interest variables w.r.t. the calibration model adopted to build a calibrated object.
get.residuals(cal.design, y, scale = c("no", "w", "d", "g"))
| cal.design | Object of class |
|---|---|
| y | Formula defining the variables of interest. |
| scale |
|
This function has been designed mainly for programmers or advanced users willing to build upon ReGenesees: typical users are not expected to feel much need of it.
The residuals of an interest variable w.r.t. the linear model defined by the auxiliary variables used for calibration play a central role in estimating the variance of Calibration Estimators (see, e.g., [Deville 1999]). Notice that if object cal.design has been generated by running a partitioned calibration task (see e.calibrate), the residuals will be correctly computed using the different estimated regression coefficients pertaining to the different domains belonging to the partition.
The mandatory argument y behaves exactly the same way as it does in function svystatTM.
The scale argument makes it possible to scale the computed residuals by multiplying them by different factors. By default scale="no", that is, unscaled residuals are returned. Value "w" returns the residuals times the calibrated weights; value "d" returns the residuals times the initial weights; finally, value "g" returns the residuals times the g-weights (i.e. the ratios between calibrated and initial weights). Notice that the semantics of argument scale are slightly modified when the input object cal.design has been obtained by a multi-step calibration procedure (see ‘Note’ below).
A matrix of residuals.
If cal.design has undergone k subsequent calibration steps (with k >= 1), the function will return the residuals computed w.r.t. the linear assisting model underlying the last (i.e. k-th) calibration step.
If k >= 2, the scale parameter will be interpreted as follows:
SCALE MEANING "no"........no scale; "w"........last calibration weight (i.e. at step k); "d"........second to last calibration weight (i.e. at step k - 1); "g"........ratio between last and second to last calibration weights.
Deville, J. C. (1999) “Variance Estimation for Complex Statistics and Estimators: Linearization and Residual Techniques”. Survey Methodology 25: 193-203.
weights to extract the weights from a design object, e.calibrate for calibrating weights and g.range to get the range of the g-weights.
################################################################## # Just some checks on the consistency of the numerical results # # obtained by ReGenesees with well known theoretical properties. # ################################################################## # Load sbs data: data(sbs) # Create a design object to be calibrated: sbsdes<-e.svydesign(data=sbs,ids=~id,strata=~strata,weights=~weight,fpc=~fpc) ################ ## Property 1 ## ################ ## If the calibration model has (implicitly or explicitly) an intercept, ## the weighted sum of residuals must be zero. # Suppose you calibrate on enterprise counts inside areas, i.e. a calibration # model WITH intercept (though implicitly): # calmodel= ~ area - 1 # First, build and fill the known totals template: pop<-pop.template(sbsdes, calmodel= ~ area - 1) pop<-fill.template(pop, universe=sbs.frame)#> #> # Coherence check between 'universe' and 'template': OK #># Then, calibrate: sbscal<-e.calibrate(sbsdes, pop) # Now, get the residuals of any variable (e.g. y and emp.num) scaled with the # initial weights: de <- get.residuals(sbscal, ~ y + emp.num, scale="d") # Lastly, compute the column sums... colSums(de)#> y emp.num #> 9.149915e-09 6.385954e-10#...which is (numerically) zero, as it must be. ################ ## Property 2 ## ################ ## If the calibration model does not have (implicitly or explicitly) an ## intercept term, the weighted sum of residuals is generally different from zero. # Suppose you calibrate on employees counts inside areas, i.e. a calibration # model WITHOUT intercept: # calmodel= ~ emp.num:area - 1 # First, build and fill the known totals template: pop<-pop.template(sbsdes, calmodel= ~ emp.num:area - 1) pop<-fill.template(pop, universe=sbs.frame)#> #> # Coherence check between 'universe' and 'template': OK #># Then, calibrate: sbscal<-e.calibrate(sbsdes, pop) # Now, get the residuals of any variable (e.g. y and region) scaled with the # initial weights: de <- get.residuals(sbscal, ~ y + region, scale="d") # Lastly, compute the column sums... colSums(de)#> y regionNorth regionCenter regionSouth #> 8247599.055 9292.902 3062.676 2933.356#...which is far from zero, as expected. ################ ## Property 3 ## ################ ## In the Taylor linearization approach, estimating the variance of a ## Calibration Estimator amounts to estimating the variance of the HT total ## of its linearized variable (i.e. the g-scaled residual), under the sampling ## design at hand. # Suppose you calibrate on the total number of employees and enterprises # inside the domains obtained by: # i) crossing nace.macro and region; # ii) crossing emp.cl and region; # First, build and fill the population totals template: pop<-pop.template(sbsdes, calmodel=~(emp.num+ent):(nace.macro+emp.cl)-1, partition=~region) pop<-fill.template(universe=sbs.frame,template=pop)#> #> # Coherence check between 'universe' and 'template': OK #># Then, calibrate: sbscal<-e.calibrate(sbsdes,pop) # Now, compute the linearized variable of the Calibration Estimator of the # total of any variable (e.g. va.imp2): z_va.imp2 <- get.residuals(sbscal, ~ va.imp2, scale="g") # Now, treat z_va.imp2 as an ordinary variable and compute the standard error # of its HT total: sbsdes<-des.addvars(sbsdes, z_va.imp2 = z_va.imp2) SE(svystatTM(sbsdes, ~z_va.imp2))#> SE.z_va.imp2 #> 1 1004846# Lastly, compute directly the standard error of the Calibration Estimator... SE(svystatTM(sbscal, ~va.imp2))#> SE.va.imp2 #> 1 1004846#...and they are identical, as it must be. # Obviously the same result would hold for domain estimates (e.g. the total # of va.imp2 for the "Agriculture" nace.macro). # Compute the linearized variable of the Calibration Estimator of the domain # total: z_va.imp2.Agr <- get.residuals(sbscal, ~ I(va.imp2*(nace.macro=="Agriculture")), scale="g") # Now, treat z_va.imp2.Agr as an ordinary variable and compute the standard # error of its HT total: sbsdes<-des.addvars(sbsdes, z_va.imp2.Agr = z_va.imp2.Agr) SE(svystatTM(sbsdes, ~z_va.imp2.Agr))#> SE.z_va.imp2.Agr #> 1 54820.87# Lastly, compute directly the standard error of the Calibration Estimator # by domain... SE(svystatTM(sbscal, ~va.imp2, ~nace.macro))#> SE.Total.va.imp2 #> Agriculture 54820.87 #> Industry 384453.64 #> Commerce 910938.64 #> Services 185327.76#...and the "Agriculture" SEs are identical, as it must be.