?by ?aggregate
On Tue, Mar 6, 2012 at 4:14 PM, Walter Anderson <wandrson01 at gmail.com>
wrote:> I needed to compute a complicated cross tabulation to show weighted means
> and standard deviations and the only method I could get that worked uses a
> series of nested for next loops. ?I know that there must be a better way to
> do so, but could use some assistance pointing the way.
>
> Here is my working, but inefficient script:
>
> library(Hmisc)
> rm(list=ls())
> load('NHTS.Rdata')
> day.wt <- day[c("HOUSEID","WTTRDFIN",
"TRIPPURP","TRVLCMIN")]
> hh.wt <-
hh[c("HOUSEID","WTHHFIN","HHSIZE","HHVEHCNT","HOMETYPE")]
> hh.wt$HHBIN <- with(hh.wt,{ cut(HHSIZE,
>
breaks=c(0,1,2,3,4,max(HHSIZE)),labels=c("1","2","3","4","5+"),ordered_result=TRUE)})
> hh.wt$VEHBIN <- with(hh.wt,{ cut(HHVEHCNT,
>
breaks=c(-1,0,1,2,max(HHVEHCNT)),labels=c("0","1","2","3+"),ordered_result=TRUE)})
> hh.wt$DUTYPE <- factor(hh.wt$HOMETYPE,
exclude=c("-7","-8","-9"))
> levels(hh.wt$DUTYPE) <-
c("1","1","2","2","2","2","2")
?# Convert home
> types to 1=SF and 2=MF
> trips <- merge(day.wt,hh.wt,by="HOUSEID")
> trips$PURPOSE <- factor(trips$TRIPPURP, exclude=c("-9"))
> trips <-
trips[c("HOUSEID","HHBIN","VEHBIN","DUTYPE","PURPOSE","WTTRDFIN",
> "WTHHFIN")]
> sink('TripRates.csv')
> cat("HHBIN,VEHBIN,DUTYPE,TRIPPURP,COUNT,MEAN,STD.DEV\n")
> for (per in levels(trips$HHBIN))
> ?for (veh in levels(trips$VEHBIN))
> ? ?for (hom in levels(trips$DUTYPE))
> ? ? ?for (pur in levels(trips$PURPOSE))
> ? ? ?{
> ? ? ? ?cat(per)
> ? ? ? ?cat(",")
> ? ? ? ?cat(veh)
> ? ? ? ?cat(",")
> ? ? ? ?cat(hom)
> ? ? ? ?cat(",")
> ? ? ? ?cat(pur)
> ? ? ? ?cat(",")
> ? ? ? ?tmp <- subset(trips,
> ? ? ? ? ? ? ? ? ? ? ?trips$HHBIN == per & trips$VEHBIN == veh &
> ? ? ? ? ? ? ? ? ? ? ? ?trips$DUTYPE == hom & trips$PURPOSE == pur,
> ? ? ? ? ? ? ? ? ? ? ?select=c("HOUSEID","WTTRDFIN",
"WTHHFIN"))
> ? ? ? ?cat(length(tmp$HOUSEID))
> ? ? ? ?cat(",")
> ? ? ? ?tmp$tr <- (tmp$WTTRDFIN/365)/tmp$WTHHFIN
> ? ? ? ?w.m <- wtd.mean(tmp$tr,weights=tmp$WTHHFIN)
> ? ? ? ?w.sd <- sqrt(wtd.var(tmp$tr, weights=tmp$WTHHFIN))
> ? ? ? ?cat(w.m)
> ? ? ? ?cat(",")
> ? ? ? ?cat(w.sd)
> ? ? ? ?cat("\n")
> ? ? ?}
> sink()
>
> ? ? ? ?[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.