From 6f97aacd06024d58ac27984e0a39d93b9e30ffa2 Mon Sep 17 00:00:00 2001 From: Roberto Impicciatore Date: Mon, 18 Apr 2011 00:00:00 +0000 Subject: [PATCH] version 1.0 --- DESCRIPTION | 17 + NAMESPACE | 2 + R/MAPLES.r | 739 ++++++++++++++++++++++++++++++++++++++++++ data/demogr.RData | Bin 0 -> 33630 bytes man/MAPLES-package.Rd | 79 +++++ man/ageprofile.Rd | 174 ++++++++++ man/demogr.Rd | 59 ++++ man/epdata.Rd | 72 ++++ man/listvar.Rd | 21 ++ man/mkdate.Rd | 36 ++ man/plotap.Rd | 109 +++++++ man/splitter.Rd | 75 +++++ man/tabm.Rd | 31 ++ man/tabx.Rd | 28 ++ 14 files changed, 1442 insertions(+) create mode 100755 DESCRIPTION create mode 100755 NAMESPACE create mode 100755 R/MAPLES.r create mode 100755 data/demogr.RData create mode 100755 man/MAPLES-package.Rd create mode 100755 man/ageprofile.Rd create mode 100755 man/demogr.Rd create mode 100755 man/epdata.Rd create mode 100755 man/listvar.Rd create mode 100755 man/mkdate.Rd create mode 100755 man/plotap.Rd create mode 100755 man/splitter.Rd create mode 100755 man/tabm.Rd create mode 100755 man/tabx.Rd diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100755 index 0000000..3ef7586 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,17 @@ +Package: MAPLES +Type: Package +Title: Smoothed age profile estimation +Version: 1.0 +Date: 2011-04-18 +Author: Roberto Impicciatore +Maintainer: Roberto Impicciatore +Description: MAPLES is a general method for the estimation of age + profiles that uses standard micro-level demographic survey + data. The aim is to estimate smoothed age profiles and relative + risks for time-fixed and time-varying covariates. +Depends: mgcv +License: GPL (>= 2) +LazyData: yes +Packaged: 2011-04-26 16:50:02 UTC; roberto +Repository: CRAN +Date/Publication: 2011-04-26 17:36:50 diff --git a/NAMESPACE b/NAMESPACE new file mode 100755 index 0000000..aab2daa --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,2 @@ +export(epdata,splitter,ageprofile,plotap,tabx,tabm,listvar,mkdate) +importFrom(mgcv,gam) diff --git a/R/MAPLES.r b/R/MAPLES.r new file mode 100755 index 0000000..63e502e --- /dev/null +++ b/R/MAPLES.r @@ -0,0 +1,739 @@ + +#_______________________________________________________________________________ +# MAPLES library +# this file contains the following objects: +# Functions: ageprofile,epdata, splitter, plotap, tabx, tabm, mkdate, listvar +# +#_______________________________________________________________________________ + + + + +#_______________________________________________________________________________ +# ageprofile estimates age profiles +#_______________________________________________________________________________ + +ageprofile<-function(formula,epdata, + tr.name="Transition",win,method="car", + agelimits=c(0,100),outfile=FALSE, + tails=c(FALSE,FALSE),subset=TRUE, + weight,sig.eff=TRUE, sig.lev=0.05) +{ + + #_______________________________________________________________________________ + # aggreg aggregates micro data into a transition-specific matrix + aggreg<-function(e,method,covk) + { + ncat<-nlevels(factor(covk)) # number of levels of covariate + categ<- levels(factor(covk)) # levels of covariate + event<-matrix(rep(0,100*ncat),ncol=100) # initialize event counter + expos<-matrix(rep(0,100*ncat),ncol=100) # initialize exposure counter + + # cohort-age rates + if (method=="car"){ + for (k in 1:ncat) { + ds<-subset(e,covk==categ[k])#Selects temporary subset + N<-dim(ds)[1] + ti<-pmax(trunc(ds$Astart),0) + tf<-pmin(trunc(ds$Astop),100) + st<-ifelse(ds$status==1 | ds$status==2,1,0) + if (N>0) { + for (j in 1:N) { # selects j-th individual + expos[k,][ti[j]:tf[j]] <-expos[k,][ti[j]:tf[j]]+ ds$weight[j] + expos[k,][tf[j]] <-expos[k,][tf[j]] -(1-(ds$Astop[j]- + tf[j]))*ds$weight[j] + expos[k,][ti[j]]<-expos[k,][ti[j]]-((ds$Astart[j]-ti[j])*ds$weight[j]) + event[k,][tf[j]]<-event[k,][tf[j]]+st[j]*ds$weight[j] + }}}} + # cohort-period rates + if (method=="cpr"){ + for (k in 1:ncat) { + ds<-subset(e,covk==categ[k]) #Selects temporary subset + N<-dim(ds)[1] + ti<-pmax(trunc(trunc(ds$Wstart)-ds$birth)+1,0) + tf<-pmin(trunc(trunc(ds$Wstop)-ds$birth)+1,100) + st<-ifelse(ds$status==1 | ds$status==2,1,0) + if (N>0) { + for (j in 1:N) { # selects j-th individual + expos[k,][ti[j]:tf[j]] <-expos[k,][ti[j]:tf[j]]+ ds$weight[j] + expos[k,][tf[j]]<-expos[k,][tf[j]]- + ((1-(ds$Wstop[j]-trunc(ds$Wstop[j])))*ds$weight[j]) + expos[k,][ti[j]]<-expos[k,][ti[j]]- + (ds$Wstart[j]-trunc(ds$Wstart[j]))*ds$weight[j] + event[k,][tf[j]]<-event[k,][tf[j]]+st[j]*ds$weight[j] + }}}} + # Cohort-age probabilities + if (method=="cap"){ + for (k in 1:ncat) { + ds<-subset(e,covk==categ[k]) #Selects temporary subset + N<-dim(ds)[1] + ti<-pmax(trunc(ds$Astart),0) + tf<-pmin(trunc(ds$Astop),100) + st<-ifelse(ds$status==1 | ds$status==2,1,0) + if (N>0) { + for (j in 1:N) { # selects j-th individual + expos[k,][ti[j]:tf[j]] <-expos[k,][ti[j]:tf[j]]+ ds$weight[j] + event[k,][tf[j]]<-event[k,][tf[j]]+st[j]*ds$weight[j] + }}}} + # Recomposes data matrix + datr<-data.frame(matrix(rep(0,400*ncat),ncol=4)) + names(datr)<-c("age","event","expos","variab") + for (i in 1:100) { + for (k in 1:ncat) { + datr[ncat*(i-1)+k,1]<-i #age + datr[ncat*(i-1)+k,2]<-event[k,i] #events + datr[ncat*(i-1)+k,3]<-expos[k,i] #exposures + datr[ncat*(i-1)+k,4]<-categ[k] # variab + } + } + return(datr) + } + + #_____________________________________________________________________________ + # tailflat Tail flattening procedure (Baseline smoothing at the hedges) + tailflat<-function(age,profile,prop,tails) + { + smoo<-profile + if (tails[1]){ + lim<- min(age[which(prop>=0.10)]) + fdist<-which(age==lim)-1 + fw<-c(rep(0,fdist)) + r<-10/fdist + cc<-exp(((fdist+1)*r)/2) + fweight<-1/(1+cc*exp(-r*seq(1:fdist))) + if (lim>min(age)) smoo[age=0.90)]) + fdist<-which(age==max(age))-which(age==lim) + fw<-c(rep(0,fdist)) + r<-10/fdist + cc<-exp(((fdist+1)*r)/2) + fweight<-1/(1+cc*exp(-r*seq(1:fdist))) + if (limlim]<-profile[age>lim]*(1-fweight) + } + return(smoo) + } + + + # CHECK Arguments ____________________________________________________________ + if (missing(formula)) stop ("Formula is missing (try with formula=~1)" + ,call.=FALSE) + if (missing(epdata)) stop("Must have an epdata argument",call.=FALSE) + if (!is.data.frame(epdata)) stop("epdata must be a data.frame",call.=FALSE) + if (!all(names(epdata)[1:7]==c("id","Tstart","Tstop","status","agestart", + "agestop","birth"))) + stop ("epdata must have: Tstart,Tstop,status,agestart,agestop,birth",call.=FALSE) + if (!missing(weight)) {if (!is.numeric(weight)) stop ("weight must be numeric",call.=FALSE)} + if (!missing(formula) & formula!=~1) { + covar <- model.frame(formula=formula, data=epdata) + factors<-names(covar) + nlev<-0 + for (k in 1:length(factors)) { + lev.labels<- + nlev<-nlev+nlevels(covar[,k]) + } + } + if (!missing(win)) {if (any(is.na(win))) stop("NA win for one or more obs.",call.=FALSE)} + if (method!="cpr" & method!="car" & method!="cap") + stop("Method must be 'cpr','car', or 'cap'",call.=FALSE) + if (agelimits[1]>=agelimits[2] | agelimits[1]<0 ) + stop("Unconsistent agelimits for one or more obs",call.=FALSE) + + # weight + if (!missing(weight)) epdata$weight<-weight + if (missing(weight)) epdata$weight<-rep(1,dim(epdata)[1]) + + # Calendar Window of observation______________________________________________ + if (!missing(win)) { + epdata$Wstart<-pmax(win[,1],epdata$Tstart,na.rm=TRUE) + epdata$Wstop<-pmin(win[,2],epdata$Tstop, na.rm=TRUE) + chk<-ifelse(epdata$Wstop=agelimits[1] & a$age<=agelimits[2] & a$expos>0.001) + age0<-a$age + event0<-a$event + expos0<-a$expos + unsmoo<-cbind(age0,round(a$event/a$expos,4)) + colnames(unsmoo)<-c("age","base") + # GAM Model + m0<-gam(round(event0)~offset(log(expos0))+ s(age0),family=poisson) + baseline0<-as.vector(exp(predict.gam(m0,type="terms",terms="s(age0)")+m0$coeff[1])) + + # Check number of events + if (sum(event0)<=50) + warning("The number of events in the window of observation is lower than 50.", + immediate.=TRUE,call.=FALSE) + # Tail flattening + if (any(tails)) baseline0<-tailflat(age0,baseline0,cumsum(event0/sum(event0)),tails) + tab<-cbind(age0,round(baseline0,7)) + colnames(tab)<-c("age","baseline") + + + # MODEL WITH COVARIATES_______________________________________________________ + if (formula!=~1) { + # Initialize objects + catcnt<-1 + anovatest<-rep(NA,length(factors)) + # initialize (empty) arrays for final output + outKNOT<-matrix(NA,0,2) + outBETA<-matrix(NA,0,4) + outEVENT<-matrix(NA,0,4) + outSE<-matrix(NA,0,4) + outPVALUE<-matrix(NA,0,4) + + # select k-th covariate + for (k in 1:length(factors)) { # k-th covariate + # Aggregate data + ak<-aggreg(epdata,method,epdata[,factors[k]]) + #select age range and rows with non-zero exposure time + ak<-subset(ak,ak$age>=agelimits[1] & ak$age<=agelimits[2] & + ak$expos>0.001) + agek<-ak$age + eventk<-ak$event + exposk<-ak$expos + evexp<-eventk/exposk # UNSMOOTHED RATES + variab<-factor(ak$variab) + ncat<-nlevels(variab) # number of levels of covariate + categ<-levels(variab) # levels of covariate + if (k==1) lev.labels<-categ + if (k>1) lev.labels<-c(lev.labels,categ) + + for (z in 1:ncat) { # z-th level + rrisk<-rep(0,length(age0)) + # defining knots and age intervals (knots are at exact ages) + prop<-cumsum(evexp[variab==categ[z]]/sum(evexp[variab==categ[z]])) + knot<-c(0,0) + knot[1]<- min(agek[variab==categ[z]][which(prop>=0.33)]) + knot[2]<- min(agek[variab==categ[z]][which(prop>=0.667)]) + # mid-points between jumps and distances between mid-points + midpt<-rep(0,3) + midpt[1]<- min(agek[variab==categ[z]][which(prop>=0.15)]) + midpt[2]<- min(agek[variab==categ[z]][which(prop>=0.50)]) + midpt[3]<- min(agek[variab==categ[z]][which(prop>=0.85)]) + ddist<-c(0,0) + ddist[1]<-which(agek[variab==categ[z]]==midpt[2])- + which(agek[variab==categ[z]]==midpt[1]) + ddist[2]<-which(agek[variab==categ[z]]==midpt[3])- + which(agek[variab==categ[z]]==midpt[2]) + + int<-as.factor(ifelse(agek=knot[1] & agek<=knot[2],2,3))) + # 'deviation coding' dummies (for each level and age subint.) + dummy<-matrix(rep(0,3*(ncat-1)*dim(ak)[1]),ncol=3*(ncat-1)) + j<-0 + if (z=midpt[1] & age0midpt[2] & age0<=midpt[3]]<- + bc[2]*(1-w2)+bc[3]*w2 + rrisk[age0>midpt[3]]<-bc[3] + rrisk<-rrisk*baseline + if (any(tails)) rrisk<-tailflat(age0,rrisk,prop,tails) + + # Preparing output + tab<-cbind(tab,rrisk) # Age profiles + outKNOT<-rbind(outKNOT,knot) # knots + outEVENT<-rbind(outEVENT,c( # events + sum(eventk[variab==categ[z] & int==1]), + sum(eventk[variab==categ[z] & int==2]), + sum(eventk[variab==categ[z] & int==3]), + sum(eventk[variab==categ[z]]))) + outBETA<-rbind(outBETA,betax) # relative risks + outSE<-rbind(outSE,se) # Standard errors + outPVALUE<-rbind(outPVALUE,pvalue) # P-values + temp<-rep(0,length(age0)) + for (i in 1:length(age0)){ + if (any(agek==age0[i] & variab==categ[z])) + {temp[i]<-evexp[agek==age0[i] & variab==categ[z]]}} + unsmoo<-cbind(unsmoo,temp) + catcnt<-catcnt+1 # Counter + } # close iteration for z-th level + + # Total column (Model without age intervals) + contrasts(variab)<-'contr.sum' + m1<-gam(round(eventk)~offset(log(exposk))+ s(agek)+variab,family=poisson) + + # changes the order of levels + last<-factor(as.character(variab), levels=c(categ[ncat],categ[1:(ncat-1)])) + contrasts(last)<-'contr.sum' + m2<-gam(round(eventk)~offset(log(exposk))+ s(agek)+last,family=poisson) + + outBETA[(catcnt-ncat):(catcnt-1),4]<- + c(exp(m1$coeff[2:ncat]),exp(m2$coeff[2])) + outSE[(catcnt-ncat):(catcnt-1),4]<- + c(summary.gam(m1)$se[2:ncat],summary.gam(m2)$se[2]) + outPVALUE[(catcnt-ncat):(catcnt-1),4]<- + c(summary.gam(m1)$p.pv[2:ncat],summary.gam(m2)$p.pv[2]) + + # ANOVA test (comparison between base model and model with covariate__ + m0k<-gam(round(eventk)~offset(log(exposk))+ s(agek),family=poisson) + test<-anova.gam(m0k,m1,test="Chisq") + anovatest[k]<-round(test$P[2],3) # ANOVA test pvalue + + } # close iteration for k-th covariate + # FINAL OUTPUT + rownames(outBETA)<-lev.labels + rownames(outEVENT)<-lev.labels + rownames(outSE)<-lev.labels + rownames(outPVALUE)<-lev.labels + rownames(outKNOT)<-lev.labels + colnames(tab)<-c("age","baseline",lev.labels) + colnames(unsmoo)<-c("age","baseline",lev.labels) + colnames(outKNOT)<-c("knot1","knot2") + colnames(outBETA)<-c("int1","int2","int3","tot") + colnames(outEVENT)<-c("int1","int2","int3","tot") + colnames(outSE)<-c("int1","int2","int3","tot") + colnames(outPVALUE)<-c("int1","int2","int3","tot") + outlist<-list(tab,unsmoo,outKNOT, round(outEVENT),round(outBETA,4), + round(outSE,4),round(outPVALUE,3),factors,anovatest,tr.name,method, + agelimits,tails) + names(outlist)<-c("profiles","unsmoothed","knot","event","rrisk","se", + "pvalue","factors","ANOVAtest","tr.name","method","agelimits","tails") + } # close non-missing formula + else + { + outlist<-list(tab,unsmoo,tr.name,method,agelimits,tails) + names(outlist)<-c("profiles","unsmoothed","tr.name","method", + "agelimits","tails") + } + # WRITE TEXT FILE + if (outfile) { + sink(paste(tr.name,".txt",sep="")) + print(outlist) + sink() + } + if (any(outEVENT<5)) + warning("Few events (<5) for one or more levels of specified factors. Check the output carefully" + ,call.=FALSE) + class(outlist)<-"ap" + return(outlist) +} + + + + +#_______________________________________________________________________________ +# epdata creates episode data structure from dates +#_______________________________________________________________________________ + +epdata<-function(start, event, lcensor, rcensor,subset=TRUE, birth,id, addvar) +{ + # CHECK arguments + if (missing(start)) stop("Must have a start argument.",call.=FALSE) + if (missing(event)) stop("Must have a event argument.",call.=FALSE) + if (!is.logical(subset)) stop("subset is not logical.",call.=FALSE) + if (missing(id)) id<- seq(1,length(start)) + if (missing(lcensor)) lcensor<- rep(NA,length(start)) + if (missing(rcensor)) rcensor<- rep(NA,length(start)) + if (length(start)!=length(event) | length(start)!=length(id)) + stop ("Arguments must have same length.",call.=FALSE) + if (!missing(addvar)){ + if (dim(as.matrix(addvar))[1]!=length(start)) + stop ("Arguments must have same length.",call.=FALSE) + addvar<-subset(addvar,subset) + } + id<-id[subset] + if (any(is.na(id))) stop("NA id for one or more obs.",call.=FALSE) + start<-start[subset] + if (!is.numeric(start)) stop ("Start variable is not numeric.",call.=FALSE) + event<-event[subset] + if (!is.numeric(event)) stop ("Event variable is not numeric.",call.=FALSE) + birth<-birth[subset] + if (any(is.na(birth))) stop("NA date of birth for one or more obs.",call.=FALSE) + lcensor<-lcensor[subset] + rcensor<-rcensor[subset] + if (any(is.na(start) & is.na(lcensor))) + stop("Both start and lcensor are NA for some obs.",call.=FALSE) + if (any(is.na(event) & is.na(rcensor))) + stop("Both event and rcensor are NA for some obs.",call.=FALSE) + + # episode (start and stop) + Tstart<-pmax(lcensor,start,na.rm=TRUE) + Tstop<-pmin(rcensor,event,na.rm=TRUE) + # STATUS variable 0:right cens; 1:event; 2:left cens; 3:interval cens. + status<- ifelse(is.na(event), + ifelse(is.na(lcensor),0,3), + ifelse(is.na(lcensor),1,2)) + # Ages + agestart<-Tstart-birth # age at the beginning of the episode + agestop<-Tstop-birth # age at the end of the episode + + # Final dataframe + e<- cbind(id,Tstart,Tstop,status,agestart,agestop,birth) + if (!missing(addvar)) e<-cbind(e,addvar) + return (e) +} + + +#_______________________________________________________________________________ +# splitter creates time varying variables through episode splitting +#_______________________________________________________________________________ + + +splitter<-function(epdata,split,tvar.name="Tvar",tvar.lev) +{ + # checks arguments____________________________________________________________ + if (missing(epdata)) stop("Must have a epdata argument.",call.=FALSE) + if (!is.data.frame(epdata)) stop("epdata must be a data.frame.",call.=FALSE) +if (!all(names(epdata)[1:7]==c("id","Tstart","Tstop","status","agestart", + "agestop","birth"))) + stop ("epdata must have: Tstart,Tstop,status,agestart,agestop,birth",call.=FALSE) + if (missing(split)) stop("Must have a split argument.",call.=FALSE) + N<-ifelse(is.vector(split),1,dim(split)[2]) + len<-ifelse(is.vector(split),length(split),dim(split)[1]) + if (dim(epdata)[1]!=len) stop ("epdata and split have different length",call.=FALSE) + + if (missing(tvar.lev)) tvar.lev<-seq(1,N+1) + if (length(tvar.lev)!=N+1) stop("Unconsistent tvar.lev",call.=FALSE) + if (N>1) { + for (k in 1:(N-1)) { + if (any(!is.na(split[,k]) & !is.na(split[,k+1]) & split[,k]>split[,k+1])) + stop("Multiple split dates must be sequentials",call.=FALSE) + if (any(is.na(split[,k]) & !is.na(split[,k+1]))) + stop("In multiple split dates a split date cannot be followed by non-NA date" + ,call.=FALSE) + } + } + # Prepare wide format (subepisodes on the same record)________________________ + if (N==1) epdata[,"tvar"]<-split + if (N>1) {for (k in 1:(N)) {epdata[,paste("tvar",k,sep="")]<-split[,k]}} + co<-dim(epdata)[2] + epdata$Ts1<-epdata$Tstart # Tstart of 1st subepisode + for (k in 1:N) { + if (N>1) splitk<-split[,k] + if (N==1) splitk<-split + + # Tstop of k-th subepisode + epdata[,paste("Tf",k,sep="")]<- + ifelse(!is.na(splitk) & splitkepdata[,paste("Ts",k,sep="")], + splitk,ifelse(is.na(splitk) | splitkepdata[,paste("Ts",k,sep="")], + ifelse(epdata$status<=1,0,3),epdata$status) + # Tstart of (k+1)-th subepisode + epdata[,paste("Ts",k+1,sep="")]<-epdata[,paste("Tf",k,sep="")] + } + epdata[,paste("Tf",N+1,sep="")]<-epdata$Tstop # Tstop of (k+1)-th subepisode + epdata[,paste("st",N+1,sep="")]<-epdata$status # Status of (k+1)-th subepisode + + # create long format: one record for each subepisode__________________________ + elong<-reshape(epdata, direction="long", varying=names(epdata)[(co+1):(co+3+3*N)], + sep="",idvar="idt") + elong<-elong[order(elong$id),] # sort cases by id + elong<-subset(elong, elong$Ts=elong$tvar,tvar.lev[2],elong[,tvar.name]) + elong$tvar<-NULL + } + if (N>1){ + for (k in 1:N){ + elong[,tvar.name]<-ifelse(!is.na(elong[,paste("tvar",k,sep="")]) & + elong$Ts>=elong[,paste("tvar",k,sep="")],tvar.lev[k+1],elong[,tvar.name]) + elong[,paste("tvar",k,sep="")]<-NULL + } + } + # arranging final dataframe___________________________________________________ + elong$time<-NULL + elong$idt<-NULL + row.names(elong)<-NULL + elong$Tstart<-elong$Ts + elong$Tstop<-elong$Tf + elong$status<-elong$st + elong$Ts<-NULL + elong$Tf<-NULL + elong$st<-NULL + elong$agestart<-elong$Tstart-elong$birth + elong$agestop<-elong$Tstop-elong$birth + elong[,tvar.name]<-factor(elong[,tvar.name]) + return(elong) +} + + +#_______________________________________________________________________________ +# plotap plots ageprofile +#_______________________________________________________________________________ + +plotap<-function(x,lev.labels, baseline=TRUE,unsmoothed=FALSE, + xlim, ylim, title) +{ + # Check arguments ______________________________________________________________ + if (missing(x)) stop ("Must have an ap argument",call.=FALSE) + if (class(x)!="ap") + stop("Must have an ap argument",call.=FALSE) + if (missing(title)) title<-x$tr.name + age<-x$profiles[,1] + base<-x$profiles[,2] + if (!missing(lev.labels)) { + levsp<-strsplit(lev.labels,split="*",fixed=TRUE) + rr<-data.frame(matrix(1,length(age),length(lev.labels))) + urr<-data.frame(matrix(1,length(age),length(lev.labels))) + for (k in 1:length(lev.labels)) { + for (j in 1: length(levsp[[k]])) { + if (any(levsp[[k]][[j]]==colnames(x$profiles))) { + rr[,k]<-rr[,k]*x$profiles[,levsp[[k]][[j]]]/base + urr[,k]<-urr[,k]*x$unsmoothed[,levsp[[k]][[j]]]/base + } + else { + stop(paste(levsp[[k]][[j]],"not found in the dataframe."),call.=FALSE) + } + } + } + } + if (missing(ylim)) {ylim<-max(base)*2} + if (missing(xlim)) {xlim<-c(min(age),max(age))} + if (x$method=="cpr") ylabel<-"Cohort-period rates" + if (x$method=="car") ylabel<-"Cohort-age rates" + if (x$method=="cap") ylabel<-"Cohort-age probs" + # Open window and plot baseline_______________________________________________ + dev.new() + plot(c(xlim[1],xlim[2]),c(0,ylim), + xlab="Age",ylab=ylabel,type="n") + title(main = title, font.main=3) + if (baseline) lines(age+0.5,base,lty=1,lwd=2) + + # plot age profiles for each factor___________________________________________ + if (!missing(lev.labels)){ + names(rr)<-lev.labels + for (k in 1:length(lev.labels)) { + lines(age+0.5,rr[,k]*base,col=k,lty=k+1) + if (unsmoothed) {points(age+0.5,urr[,k]*base, col=k)} + } + # Plot legend + posx<-0.75* xlim[2]+0.25*xlim[1] + ce<-2-((30+max(nchar(lev.labels)))/37) + legend (posx,ylim,lev.labels,lty=2:(length(lev.labels)+1),col=1:length(lev.labels), + cex=ce,y.intersp=0.8) + } + if (missing(lev.labels) & unsmoothed) points(age+0.5,x$unsmoothed[,2]) +} + + +#_______________________________________________________________ +# tabx Print univariate and bivariate frequency distribution +#_______________________________________________________________ +tabx<-function(x,y,prow=FALSE,pcol=FALSE, chisq=FALSE) +{ + # uni-variate table + if (missing(y)) + { + n<-table(x,useNA="ifany") + N<-sum(n) + perc<-round(n/N,2) + + if (any(is.na(x))) # with NA values + { + valid<-c(round(table(x)/sum(table(x)),2),"") + Tot<-c(sum(n),1,1) + tabf<-rbind(cbind(n,perc,valid),Tot) + } else # no NA values + { + Tot<-c(sum(n),sum(perc)) + tabf<-rbind(cbind(n,perc),Tot) + } + print("Univariate frequency table",quote=FALSE) + print(tabf,quote=FALSE) + } + # bivariate table + if (!missing(y)) { + tabb<-table(x,y,useNA="ifany") + N<-sum(tabb) + # 1. Counts + Tot<-margin.table(tabb,1) + tot.c<-margin.table(tabb,2) + tabf<-rbind(cbind(tabb,Tot),c(tot.c,N)) + rownames(tabf)[dim(tabb)[1]+1]<-"Tot" + print("Bivariate frequency table",quote=FALSE) + print("Counts",quote=FALSE) + print(tabf,quote=FALSE) + if (chisq) print(summary(tabb)) + if (prow) { # row percentages + n<-margin.table(tabb,1) + tab1<-prop.table(tabb,1) + Tot<-margin.table(tab1,1) + tot.c<-round(margin.table(tabb,2)/N,2) + tabf<-rbind(cbind(tab1,Tot,n),c(tot.c,sum(tot.c),N)) + rownames(tabf)[dim(tabb)[1]+1]<-"Tot" + print("Row percentages",quote=FALSE) + print(round(tabf,2),quote=FALSE) + } + if (pcol) { # column percentages + n<-margin.table(tabb,2) + tab1<-prop.table(tabb,2) + Tot<-margin.table(tab1,2) + tot.r<-margin.table(tabb,1)/N + tabf<-cbind(rbind(tab1,Tot,n),c(tot.r,sum(tot.r),N)) + colnames(tabf)[dim(tabb)[2]+1]<-"Tot" + print("Column percentages",quote=FALSE) + print(round(tabf,2),quote=FALSE) + } + } + +} + +#_______________________________________________________________ +# tabm Print multiple regression and logit model estimations +#_______________________________________________________________ +tabm<-function(mod,pvalues=TRUE,digits=3) +{ + Beta<-round(coef(mod),digits) + sign1<-summary(mod)$coef[,4] + sign2<-as.character(cut(sign1, + breaks=c(0,0.01,0.05,0.1,1), + labels=c("***","**","*",""))) + if (attributes(mod)$class[1]=="glm"){ + model<-"Logistic regression." + fit<-paste("AIC:",round(AIC(mod),2)) + expBeta<-as.character(round(exp(Beta),digits)) + if (names(coef(mod))[1]=="(Intercept)") expBeta[1]<-"--" + if (pvalues){ + tabb<-cbind(Beta,expBeta,round(sign1,digits),sign2) + colnames(tabb)<-c("Beta","exp(Beta)","sign.","") + } else { + tabb<-cbind(Beta,expBeta,sign2) + colnames(tabb)<-c("Beta","exp(Beta)","sign.") + } + } + if (attributes(mod)$class[1]=="lm"){ + model<-"Multiple regression." + fit<-paste("AIC:",round(AIC(mod),2)," Adjusted R square:", + round(summary(mod)$adj.r.squared,2)) + if (pvalues){ + tabb<-cbind(Beta,round(sign1,digits),sign2) + colnames(tabb)<-c("Beta","sign.","") + } else { + tabb<-cbind(Beta,sign2) + colnames(tabb)<-c("Beta","sign.") + } + } + print (model,quote=FALSE) + print(tabb,quote=FALSE) + print("Sign.:'***' <0.01; '**' <0.05; '*' <0.1",quote=FALSE) + print(fit,quote=FALSE) +} + + + +#_______________________________________________________________________________ +# mkdate: creates dates in continuous years or cmc +#_______________________________________________________________________________ + +# Dates of event must be (year: XXXX format; month XX format); +# option: cmc=FALSE years continuous time - cmc=TRUE date in Century month code) + +mkdate<-function(year,month,cmc=FALSE) +{ + if (!cmc) mkdate<-year-1900+round((month-0.5)/12,2) + if (cmc) mkdate<-(year-1900)*12+month + return(mkdate) +} + + +#_______________________________________________________________________________ +# listvar: # Show variables in a dataframe and the associate number of columns +#_______________________________________________________________________________ + +listvar=function(df) +{ + temp1=names(df) + temp2=seq(1,length(temp1)) + temp3=data.frame(temp1,temp2) + names(temp3)=c("VAR","COL") + return(temp3) + rm(temp1,temp2,temp3) +} + + + + \ No newline at end of file diff --git a/data/demogr.RData b/data/demogr.RData new file mode 100755 index 0000000000000000000000000000000000000000..59e518957d29e20be925696d47fc900a93c19b8e GIT binary patch literal 33630 zcmW(*2Q*yI+m0S3L=Zt(y#~>vt`a?Z5Oss-!B5XByGn#aXGIcaC3+Wa^(7KzMelWq zw(6>jW%>NSd(J#_=A5}_?(;k|@4RR30OIf7`TvI+)=KsYw#Xny)8m9lW}o>KLsqg} z0V}9qT>P{c88frUK7*h6^7?h;qmuWX57;`zcEv>F#Ba{72M5>pH&4US{&q@0*P*}J zO`F+Lu0vMhtq#iNwB;Ue-n<+tu{|9_9u?kH7+0M@dkUfEdN zGZ7+?^c!W`qgFiWm6(L%R4{pv~fgq zl7ptG1*#hdK?jj}h%dw|N)-38JY3O~$%2-ny6vCjmv=5C_oS-=KYXhw?8mP(QERAmWukce`q6^wJ&Hoj2#o zBQuC!*k}u?N28<|D%*RuJbl^=q=~1}|9~U?s=J=_Ld1mcDkDVC`tIHh5xT#FjUC9V zj=1I;NtKCwz&W+rcais#PV3f*CV`{JZWYt1XjLcp)hjh(#&bYcv^VYC<9fyaSb0#> zk@gh$yBMAqv3K9GUCl4YM)!o=n2=ldNT2ADUKf8zP`?n68hLdFl#k5qyn9VGBBF{V z2~VWI@>CCp|LW|$J`kb_+hu@1Fv*QOA|9Cms!XTe`TFn#m=%>AzFsHlW;2qnD(*-M zr9TOdNG5Fh-^XwkaRr;WZEoI&^*&WMq+q?A$nvB0=0k!40Z)kQ~*~ zeHTKCxmy$i(d_TK85D#;y?BYD_KD`c(qDbdvm`!~`f;UznNC9Ix-;R)vJS_0O-I&40e)bnwabY zyVvPM0Z}}lhycRq*N;_wh;oTk9N*hAy28@bu|-`;2Eh7e#?A75B)tC-AJKK zF8wcMhA%HzA__`*!Z$Lk#v}GIk%YZhdV84u+*jt&G;z#P6KM|w1NU! zefLQAdhbFw-F{2?x~jpQS#f~qkcW&MawG-nk7QarrLVyWah?x$X+IY7QlS;Z=qsvqFsv!+l1azr-a~mAdq?{P z!+2;jlHM=q^d-JuknZC$%rMj-P;(*gAKe1&QJhoZ4kwcRXr2p7>KErz$ZQNDT1S{9 zJ+lm93%ZZhgR5P$taWoGKYPVWd*RlnsbNSi|GU?eWvxDCT)pb{mCQoEC|lD>v(Ec@ zBiI5DEn-<|BB~HwFHSW1kr$`WZSO^0ml_{3u)yty82PC;k6UxDg|Rq5iBI@deP~ge zd0}KI`d)#`h{eH_d79?%6pb)=0M+_*CmlIEIjb*Z%=G}+ZODzH-xP4pQzBU>QwRIS zn5jm>Z7ITj-}DcfEInQyjPZRT8}oRGq`ycVnl{MgLkqU)Tggt@O}gdZ+Px&$w{n%D zM1?C-gRrg@t3l}3roV?#V0iMPRng?(Y0Es-yxx&aZC~xN7 zE0_oW?OSDR*Y0!hQN%}n|8!o<#>6&F<80p1DAwlu{> zBI$63`B$vxbZVOzpAErjzZP`(^`>}1U*DOpK=*~zBN-xmN^?6;d$4_HX%?v>*5>J==+H9LY!Fj*B$ zgLyO|LIKz;eaVi3QioZ|M1GbbPCU-25Iy(Oh)@SQD|5#BXZ~gN8(Lm{%+E%wGE`do zQ4mgjtvB)v-fuxgNF~b)1EPGr({spgq<3?)LaQ`rH|1wF`jp=5N1kE(l{?QO+9-AtKRh0aKH5sjqHW$vv7wRAW?sKe zIgZT+uBCuJW8W(cWgO+F0I4qWQ;>I?6$O#eM~Z@Un({+AMt?hBTQpg(}Xh#D{q=yz!f0o5RT!p|^ zuPSYUhMhSn^OVg_Dc$n6VEXkwCqeNa;Z54cOEF4Ao2b?rHY(*cdD!=n%RAs~bjHt= z_!!$KYm6xa4^iU@FaAfx$T|?l15BHGBL8cF|MT> zSCND6T^k*w2J75t3mR1m9E{sUd9x5RZW~~u{PbZqKNdK`*7(MhP4+(?NeF7^a)*Pf;|}UufIMFB zn=VPpU(Do$QXZXg#wsOk%cJ1=SMs}>>a-W^=DfRH!}*`0IgUJp`PTe9Z?_N+Cx*P& z+IAX(mU2XZ3(mmMgL#eW_7U1Z@jfn40e>I;nypDvTk4#-TXD#vHU9hsx_ZFL)oZ8> zDivr-p)PjckK&em8vY3?6=8;cE`Ri5yH^mJaD))lzSbo?+|%IDBQPC-xNeyC^gh^2 zo}_)xq`csU@%Us|Z7~|xq+6j$qHj=;~fg9ou?ev@?hHg)a=GC4c1jH8Td{vJ}0b+_%Wk!*@)!mI@pQBLc6DO7eVX zrfdhKXhTHJKKk*YGHT+_B)&GMXmf}Syh<4^b0;Q6Ql1;AM?~dR6x|lu6zyEG^5~JF z3i}a_4ykhH6W$cTgIR%H5i-(aM4lW z?}j#jq2FIK^l{OBIeoa^2X|@LyhN&rv8-F=r3(BNCI*1b4Pxe08S=hHm@7x9he|gD z6On}*L?^3$w>R*gQG=L_KMtDfG-R9EB30~@@@VN!l0T{B^dXVghAIBw2=3!TmmJ2l ziSH)A=Qv9Kb~pl{Be;)fbC^oq+k38JDn)_Y=v>MM6N|mKTN9F@_MvKejLYo4R+s@% zT<*LET<8zF%%Ev5=N*E@cXwlQcxs~CiZ&Qe4Arp2Cn=i~wlUY{7|`Tx))GeO#I z_;gn(I>&MZN=En(!$ntiLDN1wa+%=!@dJeTgdq>9i-`b$1W|1RXgP zlOFZKJw)z?%#}s8HR81MIIF#0vE9AcxRur>UJ227?G$&#SUp`-?Yb~!!R?xZ`?!s z*AOX&2UMOz4qlsX@XtfsKkBRYv{a=7+wLh~ygn!$HO=jf%u|5Z{2R16p`@qTQ5;9D zbJQcP5%?-K02s9Ah2a5I-}|6))HUZi0wpK7p3^$iM2G#9rCDX74TU>UTsV12?{?BP z^znE)r)V$Pi*yJ7Qe>{A36gx|AWm`})mCRgeQ8}^cjuH7P${aW+;_Pnyh}`l3or)~ zlZD&m{LkqTpq&OC=dN@t`&wP#pbv3O(AHe1*gj=IQH>|L;>L(J#Kx62Tm&<#&T-+L z!*H_jeUn-qGy@>JpiIdkZOwm}fL9iStaw6#bv^j&R*3|@&ACBe) zwIR05G|jIk*oeVv{efJMrk`SBHEyI%R5ZlNL;A9I9P={owm%O`e%=KQIAF`6*-LS- zm3l2l5)@byf7M!b2d76$*rkdJ3zAV2YfYjsE)j3#T1&|s_NE0dxV}Nea2%b?7LD|_ zYc;Z*q#}l8V>@I{QhJ9QWl8i%>tomjFg(NXJ=*XE=o`97_BBFXZ0C?5#wbE%pQyH2 z4Jvl>pLk09-J|%CeWsCaWi^}vrc$GwXar5v&VK8T@&eBrK;#lVCg;PF`ve;)!SlA0 z2F?q!vcz}Eb)wFh9Gba%zY9>?od=+8{|^&Gy<>%A#;oRYE0PC2_#PV?zx(?>k#h}C zG&{>~Un1$tE{=WqTJuN2?7I(dWmc}$e_c&_%P>@)iA+&`lC(Wrp*5m7A_}xA&3NK( zt;%u6BQN?s#D7j5x%X^YY1x`{_nYLe4wl_-x{+?oOR1YDs^Je4N~T#)t@2fDX8yFl zND`gfBfT`N!IEBb)WlwXdrWuYba$6Fa*+!|u{6J)psw6B-T(E@<+uagnWFsvP)Xk3 zAy{m8C}vT!25}c>(YZubKMe@-1;SS8TCZ!^?z9s958j<+){)?Rt6_B{tFI{A$R@AVTLCpI z<{<$`+FVdZg5o{pg`IVBHLSv#Z*^uUMMx9&Z}{Ez`+c!_xzi)@5m(H8e^kZD{Y?8- zIG$S4g5`GyB=&`GeSd06p{Ou3>R@{in!cbu^C#SK-)Dv(Ix|=eV*>}&s3XQ5gA za$?S#2gm^s5d<+S_6op? ztiGDVx8iuzkXR_|Y+X8%RRB-0GnLaS;p-M(=`O zzejhQdHm{o-(4P@QenS6=V{{eT*3XbW#6C}QfTYH7z&F~!1EV-PY0hIy0huA&j#$Z zwe!J#9={o6HV`r0URSddv%JtP{W#QKO(!j|)fnkty-8b0lWfJ%B3LQ7+F1s}Y};vq zFV^LX%h*5hgFQByIK;Vi6$(3x98QSS_N#CTi#I~L;=yGyJ)z3F`LLExp&hW%&Nbnm z_W4t+=@msD9}Ztoq|Le;W%`^H&~#ZY?!x{#KjAWWz~tCaP#Ej^ZB*)JN{1N1(JE)8*4w(frUVk(jn@UzJCvLcXcZ^RSfQb@=_%#}Eh_>H}+gKJN2 zt5^*#mQwXgK3$5eTGxwxW#>#dXMzY|nLK7eZzBJM@N44NW`WW(*`v-Ce1S;z)F0@L zA4<0J2DR@0NhiJpXn4%H1$bF4RB1HIoNNXX^g8?l^;uldL#>j11yLL5j&8ySndU|F zexs#wAm3i_xcAFmF4a1xku~4AC4qJ`(st8n1>l zHw%n~7;g{dw$}MfF>Xl%u_RNAFj&xzal2x5;D?iX)-wfJu%3G+-avfer)oi(yh<*3fVumR;*7`3_Uhu~g&a$z8{#r1`(y@D1!0a?g&KDM@!>oy~Ud#c7v@cIu z)yw0uDs^X}py?e?1|goFZOc#Sccx!$Z6q(2h)KSk&T`E6KawBuv!9HF`GySP?MJ&kkJ)ftxwW zDc_0^ATEIFb|=^$`4+C;Z(Z&p=Og|DW&z6m(DV4dNyz@$v7ZZ@E-Tp6Bld=Bz)cx~ zUlC|TkQAAPTn$epCyeQIa9blB@?wa}{5arcd$#FmVGrN%+P>4&^~_K^HZ5iw!Q)Qc z8`D`fDuu4cQfsdROMb}A8`M&riCNSenZaMzSM)1!n1o;Ui8eUR zc@G%yP7OQ|6o7P{vM=!WsPf3um(9t!+WUunB z>QtkMuZN##9XI(|C&4cBZR^gwnd4)(CFv&~(ZA0jSXVoNmr^5_&2`f$5kom!o6Kmj zVv=z_E%SNh@~uf@$BBW$?zRQQ%Ll&stu52i_{CjAia@RV>#c3b9nts#&^@H|yyO5p zZ6aRyMR;3WTGmi;dw;qRf5UhNo?h|!i>0Mmjhz|Q=d(kZ-oyI_0<$}KbJ~~Uj9W0~ zc|eJ0eI86ncWm9Ni<^Di;oZm0A?{CKh1;QTCSInAmRqg<`ApB-I_omN5L!;N?P~qg zdp>R=#7p;(w+OXt7D6)Q*-vA(^GIdSsgr8lzsbRbBN2nb_Z_CC{rq(FQbF$N zbb?04J!kQ#YxrMHg(}R#)$jiNt=PTYs@WD{*uv%7oC4fQ5c~siyy(G0gLyw-x6s}g z^!%uR3~0SZTP5vl6ZzI+KM{>d&29KW{v~T#G6cliI}~T!tI1n!ULz zbCrK;S@V&V0(98h#d2lbs*6kp_zc-9lyoCOLW)0*g4apy~4jTRVfX=hcFi?J=iX-I@#i5a;XZN z4~ltDH3+x)YX8L`Mt`Gnre50OH~0GQ+)sy=d_9Zhp6Q`iut%BjsUMf?5Air&raj4h zaXK3<)L_O&7JgD-T_op4hMnry@qB^PoD{9EXR8!pD(R|6x{4xUOZJhq36U%8E|vQEv+dDwe+YIn zV3vU}QkUI&=2cuhq9yC&gYK#oigl%V0pd$#U_i66kC`eeTfeAJmt)=`ZxNm zT|?*a-$s5xMOlF~l}F8BkdQOYvpLa{pkUNYO^cfmF5GkwKH=m(7GSE(+CCV3$fMEn zG>?#IX`thEnZJ`IeD-{zrFmTRa6WYGmAMj= zaz3p)Kj>P%wD;whMRkz-+rp!QhP$oHR|k??E+6h;a#uA_Q)i!$VGGa87_o69?cutC z&jc==y7ba4Z>rCbD1TlA;Mr;SIk;fwm3w{v7NUU7k5@tNUz>yc?iNBsdD&|zF}IxN zLe`82mj%c_NfaOUKCjpMU2oRzl+R3C{5Kb*HNAh3HsP0EN}+C%X1cLZ>z6)uY}!h@ zTJ?BO)gy4+ruSh7|6o}HD)Q&Yg65gvsnFh1*fHwjmtlGUswR8Yvsl4^#XdQqF-KmN zMJ9}nsj-r4yAb|jB#2gCXR9=aWYy8;R7vMlJj3R*pG@Rsr=C0vYt72qwmbt;9r9I5+PS(XX|T}&PFYkS9)c$y0Il%pxea(;sf1z|@Bd5W}THyNF zlu1u<)5;%o>&Wnl?59gDxzdu*Gw8&H3wLc!LI6U74PV7)o?Oqs{~4ee`f12zu|csa z-4Ue=x6Cv&aS7bfklML4+xSl$dsb&1A7DMa;sj6r`@R0+&pz#0IqOK7qW*M|?OiT& z-Ej#Y;TOy%mu)pc0fm;Rag!v$Gn znW;~$!%$m0JYHTa$y(CNS|@X3y#>MXNU!Z+0WfpOtXB2M4j6cIXlp7xUa&IsvQC42 zOxGEJd(wkvD-CKQ^RkkOFsqqt%A0Jz>$J`nisI|3>~#0jR|GIyCY>zQUo+?*p?9g{#ZTK2ehP1v1neG+XSG zanil!t#Y6B3_bvU(a&_qwJcis)SAdrdedVN~i`~S!kMd1L z80M-Bjyv}1>Zt?-5r(t3xN|ZN&aK~eBMN$tG9TO z|7jIxd6{VKe20WK<;c?uyje{j7ed$`swD^S!)BfBtN(QG3Hi^S4UYfQPA$?`ml2tz z9;0Dta2PgYNSz9|=BN#iNvciMc&XTOnwOwb48<+r;`}~<7N-ve-=b#42ti_Cn)O6^ zKlEPsmXgAG`CKy{&+`efw@Ao7p{PK)aL#>d{EPJdparBtVlF})VeKPS%3vKbv^ZK~ z@HSr3nl2Eib)h|paT|eL3`k@KQ+RWObzddw%IHtW6l8z2>gj^n1&?K}cz4sxr~YRb zIQAr}-rNsiUlI@v8|I&TRC!vQFw9IYEq!fJ#?J)I{ldmO<_-+Pwmejd9kVG2tbOY{ z7cuns=TDD6_3Eg$XAYnpE)h+I=8b3QD_%Lh#ujskVg}D284+Y-BZUxOL8{jPAX#X2 z@B{Antb)h95oZy4SdVK{yEvl_RY_irVd@B$JMP~ zsV|#e1@DszP#0Ml)gRkjU+=YvXYy(dD%D4m7+@^lobfd-+G3`wG=eOKWy zdz%7kDm+Dl!a~}EvI$d&jhi8*#<$8BLB|5KhjZgm&P%-UM}SxB^Tm>Ht_hyupH|=S z4t&Gd^}8j-u@sfNe!h8?_=X)7W7*%TTRV!Q&)QUoujdZio)mCuiTlDR-Ex0+o+;J!RVZ{o{(*`PeZnNN_qoYqU+WMfb9#nkP7AA}Dp> zoND+vO!=?d%0<3Je?VqqvuK!JZFjk)!tS1WZSp5EQF}Y?tu?Dm73o^et)Yz1@CtI{ z9^3NuU=VF#f}vs+#O=P#hvl8*arCx%OOp5jrA3zmqygASbGh8d=+Z0-=b7pw-U)WSwyH%dp18;xA3%KC6oFKPp}Dm&snpg#oBXP_-FXh=c)DTT{rKN zv*AQPCLl_8+9rd5@q(!}%Lc|HSHFZ8lq<^ZuQ8@0=SG7nZLl6~b#YKO?U=BYG12MJ z)Y?J)RhiZoniCeU>IEuabGZ+Gle_j_iZTr!JK8A8GT#2&vi^oBT_>=yAIlE`>J-lH zS^fE*t+NRUF?LuFb48W=N$vQyK1jW`QI-$ae)|Zt`~uh;i+H+MC^lK{cZu|Nhk@Rf zR?WDvXPgJNq$-ZX{ifuS+K`t@>kH7QCkc!HB(JDJXPn=m9R!9kA!l4VIX*t9IX5bK zm>>6cGJ?v}Xp&5>l<#xqi&iw(RybI=nR($aFyJB6?o-o=4jVYtFa46JUg*n!-SuCz z>_a9dsFgTGBH<72W6-cy68zckR4UWl>wT7m`1-*t3^eD`X(k0pE~=)YDyM4@0YE~&#HoOgFI^+XA>i2N^BI|~PgtB}c-kZk8c05KDzatua z1b?LaC;4(6Ahr0!adgYxW^?$Jxh(BwhY`#6rk=dSipWlw&xPg9eZsxf3iS0a%?1eV z+{<(-Y+~uC1K}C{X>al#lO^!yzineDVchJo`2+z)>ZawOPme+xg+KQ3SY*T}4C zeZrRfy_MG=5fsXlDYE|Y5Zk@eE9aOqLlnMH&qL)yH zo!2$Yq?rN}0?nb>_pZ*x(RN7XKOYN%R6?|V|GL3_e%A8G!~s{J$S3?5KsyzL?R5Gqq%%H2qTrgnk3^+@-gvbRcX1TfBKks~eDSlme zL!g=3JV=;5_Nac)Tk1AmI*WJ=Ffn;=o=W)jBJ~iyb*Ox6KAEaI5h=g7 zaW{1%D1PS`e5V_92LGkontVO)$$XorZm#_V_T8C9&)m>-NL|J(g-v3a&tm7T&Ru5Z zdb>lnTqjqBk|$=J`+TF8&L+;}auG_&H%L4hjr-WP*0INr4J^8H+1{w)8y5VkZbnrW+TRB?T^O51cP zvVvVUw(j}d>-IaM`}26hO85(C`}M$;$P@V$d*B&Ue-tA_=1b-0wZApZd$W1yWin~4 zNG+T6Yl(Yhi-5uD3%hSk=6$rd%W$F`=hEdurT;*gK2}vfb&6%Hmzpm3NcB#iodxO< zgkehMb}!4PsPt46DmVVcCBxjmD&Rb;H(<@~8)R_$MsAm81f8_U)3jPcF?pft>G}k6 z7f9$r%k=3SY*K9Qcs22E8bR>Obj4U;H=fpfe!hK>jlLOx_d_K$lqUe>V_{~(e7X3y zJU`SzPq)qGjw{ps&*V#sKeO?iRf}^^cbC}9LjBXu*iDgTV9)V0Z|*Yn{o%G#Ja=ii zSZlk);y2hh2vSpdnh47QJ(-t6u=CkHdhoerlb$bJzN7d2_cPN3=#pbYMigT&TxYw; zn`d@X1yLn;#%YC$?dR~^HrDEvZy(Z@Elx1m4y4|gNe$VQ$Sj1^>9B8Za@m-6B#?Uu z5L{aqJ4oDvTPX-NKN+{yzF>htXi?>4%e{uO+0uvHWtMv>9I~FeW8T4tx%O+dF^A7a zHtTj0d3A4bik9yzLzWg8rtl45itZqTS374lTmPAY^h1KGOW#j?wSif=tP9p5r>|V2 zU-H&*D72@?rj31O&s9{=a-Q*l?Ib5d6i-eW*==8J=7QBeQOYSBEqZGTX&atq^Mx7+ zt!(8NbC`cDDOra=2a22v{+6o?W=HNvSK35=`o}ZHyM}AF)gn3F6R1aq8#FWk z`JPnlhi&4Y>XkFU04KYed^1W+XPA%~aO>eeG!w05HlEXI)yg;A5p*az^*_^`*ImNw z_fGoirG*-PTDlC7aL1qJnjLN^kG;S{l*64fBv}P_>U;K`F)nF|;r`#7@?&*G91%U$ zdu>^;dw^cI9a-^i7i^I8yGDWaxZ`qbBNf##?JRBKm+`vvtgpvP2G~qtPbR-5xvQdt z*R#UTGgS3rMGVoNsC)QgG?Wlj{3PgYPV0QkS~LRZUON@FmIb2q&8{^rU}~sd^!moA zNtKDgXm-C5iX^387rkkkHqiY zmmok+V3uBZ(~!gW4MxfJ`4gKl=B}Chz)~`TYX`f@U}x3QvSj*7Lu=a6d=j%k=#z); z;mrjXhzDVda>ir*nl5J#H;nb1y}p!r3H)5tr5;y&KUl>(_K~GBo2uu8(>*>gO&{`* zSNUkThOd>Ni$9t{+$7-A_3DuC{eEin%~ ze~a#oZ*cJ?Exj6Wnwgh;qhcKL<8Y5~t-@;;;@f+T=g|q__gblIFTB({?h3?C`Mh~m zZQq*4o`tu8wrdL+X<2RP+MD8^YR`iMgUs-w?*@v153bLy*BzGg7hyGEtTu>s6N^r_|4dh|s!=bKx`uhqBAm+7fxmjBVAxo_XQU?|2j7*bB(@P0_e> z*FA*D;K0F2y%#$u>l0!4(dWRK5KUI+?6-M)d|Y~;u4-CNQDJ)chK;cn`a|+BOGx|# z-+Wfu$#}|-`BKs6Z#s4tD1!$+E9;%z$P%%kp}d25NL z=&+e;6_U};n~w{LAl_XTic0B>7d;!3aRHP*cu@B8^N>KG_9>tBQm-n=x9_hg*GN5$ z_Us&;4nT8zT^T;b^4fwG7QpAyhLh3P^aZt zb3s<>Lb^qwM;zqFj=;L=>(s()Hq4w#;n?kySa$0bfuj}1HQ!M!-sCXCN zd&o~u5M7gpf)-wS1{Cv`%QB`Ofv{0YbvBP0wUa-K_km=B&pgsV^MU`S&f*I%yGxsV zKmfB(;mMdoxoFL~PVd^&W*=q^!OFmPyawyr2kxmi-qo62HPJA8g0HdRYG3W)%Awom zls3^b+QO!O<_{B858TBv$;FT^Az~2)!DHE*hQb0zGVq~u2ra8-b5=WZftQkMvk`(~ zS5_p^rVt{vlfbTj;UDV*cB)UQpLFHFu4qj-3Z`Liy|Up!omm6JZemvpwkLiT2nKWr z7K4Va{(>^f@O#E48#x*5t71z6iPHRKvZ72|*+lxGvNl$h6|p_agC8ftS{ih0^B)6HmAY$dBblbAQauz8MIV~@0^P&b6)AL58IG z5akQko>$lntPacdv5VG?DX{etd!9a>Y;4i{6QBiKc#wL5$6L90iIC=ofaZKoPeWyV zQbH_w`j7Vx3w^yeK8WF*1P2t#5V}o3XD{)%dk{{Z(?3@1FtIG%oI2dVM@}<6I-0o~ zv6{Fspc0s)$@#}CqoAeC;mN2e<0Nc%n%`3A3lEj0J8+fZ1y*%E0xb}CC0?xD z${(y6Y^p=|)D`6q+hkujH@^u$&ErBfZr+KcPGz(y9DvGy7RVPot81_%a<{K|g!!-J zLtiGm@!8bcI7eFVZ1|_yfN%AyZPq%~GMsy|U(4IKn;hiNp9m!X4az1!l$~!FZc6$! z<+nP>_WN^iqn%0noy~VjO}xfpz_Bu(Ga1tU{4WU`%(fq@qZ`q50}>_Trx$0A@g*jO zhsNosYQ*0M^Og93)8>8pp zmrvmGH@=(oa?yfPYoQY#5|X%!WQOr!!AAk6{Ph`Lx(b4030rV)K<490RhD~G7HnD@ z+#W?V)|bLRMV;xS=&BZdA2|D5G`f2>q{Ey!%;tT=?mK&-&Yy1B+of(md_!Rk**5rw zIU_xU-?${ljR<$9EA)~14ASSx`lI0M!!`I(!*gm|IutoZ3ctDWa#&th|!_R8Y|!d?w}nW-e?y>IcaRQ9Wi z3;6rG7lLhP5Oi@J9-f+~eY89+g7mNcYeX-><$~|Nj4vC@40UbUYyO$hQEC0T^D<#7g#S8B!9iISV+5X|<4()TZqv!qX?+{IqB z$oxa{7xdXN%{h(^jWyMEWG((yEcA6FAnZ6!=_Lx^yqEU+_i}9JtYj@?v@PJeY+bu( zkolBFiVxr;x`S2ws_4=jlr?gjdKb}QWZwnVjXiT)oaA z7L!ItgXNtOW*DLd?jMVD6~tP36IUiJ^8xe= zmUW{)?Ozl!Cc2k!wDzVs6Eb`m`wAGHsS>Nz=9Yfg)Y@g{_vy@cMo=aWyPg&WF6Wxk^ykN9=sFpEx)bzOi&&>O1(>3j5Zu|>R9|NP(?HK<1Sc`OYyQE5 zDyt^AVhS3Z9o>~bS=fYc0Bm)A<*1~lWScd8FSKmBU1bE7lcjge_SA(QB7v_$S6osp zUg(#agZH1CW{*GdnNbxP+wSr(8~heipfdJ$d345=o1cdED#Tk&^b>mNbfgdr6jW(> zw*d8|O~2X)HfAL+@`FD|=O3+_kXR#c@9E0~x8;U>W3+5d{QBT&^w<$UV|aXELBPPj zzMseA_b|!v*7a0$H!lE2_-zi$CFh9RG=-SN{cwTi3)DcmOB)xVs3OpAkTbzacrg5*K<$sIsrT`wirRngr`XS%Dz%QXIww;Jy}es|oJ zTe!+v`EK^2aHbHxwpsKq(Nj@p}WG)(dm~3NAk#y+hJo(akZ-WBph} z*(ce*V|OA(dVaN5E4yH$ zbHMl%tVdgfHK>l8L_xVcL<7+}uRHHFKVu}3tfQ#l7JOqSMY9%Q5{P22IQ>*Bl5kzm zEAe@({HtA`MF<;et3iK0tQX?^PGn3$)Si86wO&mbzc(#D{lq4uTe{6=3xAIJG|6UZ z%4`>A9y-czihfm0HN&7k_mGEGsR68EnR!Q!47DFM0kf;SQd$Z^=?{?~4_UZ$I{`*&iz0 zUW%(g-8T7I&8BL?M-^LTHqhPrwsQ^dmcRZq`RB*$K5~i|;Z->fPTKx*ss7KjoG}** zYto{^D2FV~V$Cc0I1Y(J_UPAmaZ%TK$2B|*Bi_Y%uGIAxNDoL_q%B&6Xfe-+dhQ`U z$t9|7=S?JTit_a09p%jYEFSvsPl(l$R13`GQ1#!WSqj^N5hhb{r1A$RmZUingSdCj z;qe`sTU|X3%3TxmWg`kth#e+T8>S7Ts2}!}S{nsU9fu?xK381i!p^g zXthdk3~pvcQQ;SSvDtoO^%DFm*L4SYDRU=?qV8saan6AN?f?9*tyWEu-8fnKYMPR8 z`~i($_!a!?%r``^58R1uAD~SH6E+L48xt0dKo zBg&q3$bNwm)aatv#lt{29($mYgil-K&?9Ecgd5e9Oe9j<@-M~tx@!GAE#a8vmjV&k zs4vvP-XEDQ0a*&!sZidi+$60}pO1ym{SM~n+Y!Ti#$(rsk!g9#RlD#!b_f-I=nKW? z1ATjVNYwsEtYdVb%hPS&-r$B}UWVULX?q#{U7(egn*9n`P=O0BL+~$HW9Zd#y5~F+{uP_wKY9KC|44TlG~%dP&gclBwT( zX>B*QD2Q69F z>>mdQrV!+@ZKEAsyU>zeqE#M<_+JKlindWYwsYlPHv0o@`^bc;zjEF+sHyAi|87fL+IonTDk_k8s8UZvr5X+)sYRuh z3hmD-Dv*d&Q4wN5NCG)jsR9CtXFwpSr&8q*ku!uuj!}d}Pzk57ArZD1U=w1FJ3IOF zzwZ~%{O0++x@Vr5d-}{|zVqFm$%{Sfy4IStuD#Y`gW1}%=ObJyrLOAHY*Wl-L_p1V zmCq|#jXp^8|5M=teuCN^O*uzJQ1`zk6!e|=Eh z81W+Kh}+Z(;!@pZ1@PgyzL0+5O56uiRnpo+Pk3EU#&UsS@0SBo_?CzGu9py#OCC-S zWd1x%>lWBRs#Wg2yCATS*v3fpe&ql>bFjTesb8WrtD$|vkGwumx9k3p+}+mRbW?r8 zbl>Z;Vi`BYdm$#BUCv=dVwU}$^uvhI4LE~ur(VnF3o0Z&NBoZvjqUZV+OiGSVOMP# zZ^5`rN1>1B6LG87QCOV10$jkU`6SlSp9|2#8_O}Miv&&xxdcj@vW1gT&if@=eTP7BiCo?=A9PE_mk44LtN`J?BHFrBcbLoa8xo{`)py7 zPg!+){n1RQMe{eyxsmB8J-IZ^)<rv4 zSd(WjKDefQG;to>=wXQQD|smj@r^7f-6Hwi);Tr1(<+1ake~Y8E-Z5kJ<1tpo^x&o z9Mw(Q+zsy%Gp$}BEw*(Q=hf}~%QPFM&kF1+yfY+VijC(a11*@&5uAvdax_A}v%zQh z6q%eW6!B*2(v2u?rei)b0TFT)*1|QH)!H!~9`ZB1NT@Y~R*l}(@BUoC zJI8m!t_Ipc&G5mTjpg>)`lk#h*|p>-65Y+B8O;m?B46>SaUA;s9jF2Sl)+wqPm%@n%eZ5pf4xD0Cu zTETOXqL_TQR_gW>vSOWhQ<-q7M@rwQE^m`HL;hOisWaz8^YB}O`67E1)7$=>@ln`& zpinbP$LP-M-MvXW8t9Ee5GOmqDc z@N0zrhY?#Q{j?8bYlqkHoTWekb~N#lu5LRX7LMvhkC-wWNVaQgoEW=E_hFwA-;E{+sD2R) z_vD#q>%_mba#pT9q5xNsCir19X3Y@MQ%bKrVl`PiwE5d)v@hSM8XX9Bow<8C>;qsQ z5##ZlNoK-%HHxeIwAmB>M9~8NEDL2svd^GzN-r2w0;fYJ1uz~WtLoC#xGWcNS3jps zS2av%1q8a&B8u?|llmdB*dp?r^?eHdfP3{Eo-bH!f(QBb1a@+i8|!CDuDmt5#mqEN z3^CYA36mmKH^EB;8t;(NOz#|pe1{wAZ#$lsZmgTcGUXq#LJD#|v$dEPhy|+Jo6v{m zp>NI=&^_DPt?*hfV9H1Kg#Dbl5lmy5#JjL$uYwk2X7I?j9|aa%K2E`3gXe6 zVTg{NCA6!bX@qp_AIJqGRxEK*jP?=L1eqiO9V`;I6I#rPTLfIBXC5+h5tG4K&n^BF zI#+vB-lanPdQrnh$0(-R`u20bg0wN619;AxtHy6vNHv*WQppTwJjlf4CS@3J%{UK- z_>jh#s|5j~EpijS3mHjq<|~-k4%k7mdZvx3VqQS?o+R+J)|yf88;X6VVmfp_DQGj9 zZcqGl>b%V|7gWD=4OWcP=&WGMAp#Iz@D!wpabE!&{zOCsD*qv>H;hetLcL7qWkSth z@G#AU=Ig%9!UNvKHki9W>dX`euomJymj=oL(v7Twoc81<%{>$PW9kZHzBES>;Ae>9 z`k8F*)cd)T?;CvN9_{{*H0QM{v&`5mnoc0SxMp10=Ynke%-rq7kzvIcQ0F0crjf`_ ziZcb@a2+QrI6BrAY5$NX$HDyBbV>`yJk$aD$>iGPXWCm{JjdCDAqp4v4*6~=hk!aW^aO5P z475Dr9Cb4dUzQ(YVpQb}EK#_3rVNALB!k=g zXm+7A)rgV>XyN%>Y%6`ap)FyH!HFg#YFS@cdGc#V#R0l#Q3o3<26r%9w2K5hGb;mS zmARVN7;4TfMTs5NqOS!wA0FkLdAwb>7*LFk3gBwT@dW~B({M))L>VQ93mrIteAHWy zyL3L$4zgP{TP3b)>k@Fvgubl6qKWfn82m@bxwj~QM?XS1WBXo2Ibj5CZQ$p8`%(7C z@B(-fc&yFneZm9uXYi{yWuW8V0!U~sbOMigb3yAGK}#)WQehQeglgK1v`kw;NEy?1 z6Pr#37iD~0>L=`9G%z(`7lU#Nn96i^eF9D5jVru4mwOiR(_`L zJrj_E2jbHDFs416tjxhTL?j0)ri&MjGg{|{?2k$|SkYYF%>v`N=`yfRYg1kWm#H?Q zCa)c3qBu){`yb>1D6NkKJHouNOp?}a()q?-grYK_>5I^c;SFc-&$VA!+VsAWwy+gt zThjc8iDu)c1)-y~yZ8`yZ{T)w|5yZW0MvvDOpcrFX4##5q6@;1O}6yIjsXvhRgf0r9?7Dz^|tEGU5oNnIt zo(FB1ivYN2nBpU8)i_G7sUR8jq)p`1!LUaX#}c**Dq5{ro~!lMC+aTkEf3WmSwPVH zag1d3m@&}Jo(@tci$$IZ4Sk+KK%SMSdGMwBpsZyEjLttlMQY;m!NVdDV1~%2uw)v@l?gLYdXPt zv50$DQlOY{oV>&;Fix=EP?Vb>+x8Mlx$Vz&1CFp0#2LUR@7sO5V2?pH{er{+ zY%>wEpom(+Fkl)!FJI;sYTPX`xwS}mqlxF#`{rf24vqqaXZ)?xImWs1X+&4I+N6U1 z9Lcc;w_i|GZSJn1<08}y(@7p-WoH6aD?>Ja0Pu7fR?oSmD_=0YG=yyj+@ZU4%n|_F z!IM8PNE>C{GHxR9&hQH$tasdrPWC*w00Yvbc%hG&I+l15*st~JCdnm_NVAMoAYW_F zpL)j`f8=QlK8x4GivH8qLL^Oir=u9gVKkWf744# hGZCt zvcO9cwX4A+Sxd+^cS<;KYEyv>tH3b5_Eq&8JnGS_(w)>N-t@0OTDfR5y8py$PTptS zuz%T=;u^V+(fUQ^i<<-0-#y+N#P$FD!+@EmOTXQ^Ecwor9rWKmy0tVbu=(;YpV5y5 zym;!&!^ii;vM7y^rMC=8DEoBTWYgH7v4E9S8d3c2!QuT?r$!08!_%ziS(whiBlyoW z=e9_)fxDRNlW>%c<;~NE`z_3ZZ}Uf)w^?ZQsRN!;e9e6xm25>Dw`EZ`nl7=Q$c2z` zsrfhkS5toNxNbpHatV1WO=|JmAbje0z33EFe2^Cdx}qbIW#fCK$QK62TgqMXeESx% z_uN$qd$U~sv~TTgi&W;lhOvE&Ci?qY-m4DxX}hZi^@!Bcv2tAfhz<3=TU<+6v`Hvt?;O+28_xW50_4cJ6yZzy{xQ$3}5xB9{AyDAQXDV98JW)LVMTxID!idtQ_vw z5(Ws1Y15+NqY{E&NII<=5iKI4apO2_B=0F!JIyTPy0Or>f7f zRpl@otu54p<@kPrOFC9Db}9^Y8puk6%H(jLu}5`x%)zil6b~+d!^4KqikcEWD+Zy_ z*%#WhI^1cyGt<*>Z3lC((LOwE>@nme*AXDROhUe4FDGjIT)xyMXl89ECSHucCD{yiH z&fs6oN}e7WuZoer;MVK>431}~1~(b@R;D8|C#kMm{7zQZ1YTsoa3+*RGs zD@hgf<-9PIj=?2XNqEh?p-?MU0Z9|GXNF3_vf&XkH(g!l+DDsTHfRf@CGi(&3$AOR zsW{AOz}pfvhVHJ5*B(KZXXaac4V!9AN5&R+zpl_4`>@Np^@Hs?r;R;}6KqZkN4(Fl zNA+n9@Zy?YPGGCJpqipyD_LlsdiShNEht;fe|fxH@!UyC?RM~+c7=Sou1*T)YlbEA zfO9#kZEd~;o6`5FBxvVkO-4HM&8E1XHr@7e2(KzE#)>)V2wW}^`s1?AX~-v=hMS`C zP)|qO4P8@}`Yz`jE62%$pi4pY_2@dA>@AIGYxX^g+0)=TRDT>+)q@Yqi6@awJMqAY zorA=es%qlotdsd7>j9?gonM%;CE7bnwAk&b#wde>WWg%}x5^nV}|)Uhihy zFkC2b!SyZ#YBdLN=r+xgS5ZOG{ScHoOkOI&0}+BpJaKVozHyY{l5VS?Cl7i;Z!Lw< zsNUf$B&D2(2FaNxotQa+~3oQdO!7sw=sNsOX4PLqWXhFN# zmrv1H7ONTIF`y?J@dl@<)f3*WXT`1U3!Ju0r)pjW>b`K2r_x2FKpZBXYR73#C|)`~ zAFYimH!Ogw)phrx+?K%KwH%tf8>cmMXk#xUAW{gP3Uctc`Q1F6iqW3KwjN_ec=3f^ zDH4p4RjK=(4wgNquq}VH*9aXZlgE{kAr5+Dp_4RW@uoXW6&!_E%pE5=B5tp{2)eD5 zi5=S$7eV25Z94BbrDmejjy>^hDiMU%FrY({0DN4L*+!%%a63+5N4b}S2 z)zjZ(3V1C#(#GKD@IBj07nq+O$BYRJS*ydwS$`jM)y3M|8UbdjjVqXM$bI9yWv?iVU!5fX=_!oEbIj-_`wMpt1z5HxbnI@!`S`%xR z;vduh2F_9`>P_isGS_k1471W@+O%1oW;%gf4y1_%dAwW&vscH8 zW=VI4v)Sr%iJ#HS)~^~DOj8MKr-r?Gf6L26dtug6VN2~(OLFAgJ2=TDN$O|$&}iAd zz!bLnj$bf8a6?-9bbZDYh&jaRTg88Q(VLH&F4F*ytbHVF0E&Nbv*oujUB``(_eQI?D}Udm6KZ%h>9Ef(QBC|sI^2>d zs?d!4#o*Ak$i+2NnY%#yUmR`c>?hQES>$6+sQsP3>mEsMzMjs;A9);{fsk?NRPAo4 z^R`WXn(rY0us(5Yu^OFDJ3!om+I+{45*HqWmoDt)c=qo4X`n6E-@Z&kK8QR?Xj#-P zVLBs=DjmuM=k!a`>9VP8SDs$u6w%mNi@Z*3#Xs70)A>23H4Rg;ns5IDTpy+UczhXR z8iSMY=?6V7zE$iU@kvU*79U)sKn_>?C0WzSkHDZh3;m3e?}{Cjm*RQ)?FXh8;dEJ= zcrU&n`bu2AwlXMG5@;WapbX`$PLC_Nq<@PFohqaltEW5%q7}B=n6ngw)f$);5{Z8H z4hW(jM*??p#6TD?FgrQ^gpwZ?DB9$agdC^pQnGS*`Xx{808l7dilzJHQO3nb1Z^OBo^;`H9 zb2j+Na9G{WR^R@VtT5X|FlJnlv6P(mMgjrKu#|0g7+G(@(8jRsq;YANCwddc$piWA zrYnxk1;2?>w%os;Fww(^Gf&8Iy|K5gBX!N`CL5hSTF8=MEr|h1Zj{4zbx;|gzi3{F z)iMLZQG_rag)|f;ej-Mx?)#kgZX7lVH0<9a9J9Re*RCJv$rx zg)m^;)tV5(;yB3&1F z@s4tmv)E>^ffxNa_|n8td|S3Us<4zSn2-i+L7jQa&Wd@@>mIg0^S#4jjK9+~5aS%_ zqVcBZN_Bm{?b9;s8T*Udvo(1hkK@ni+HDq#%08LY17I(J^Y1E?G<0jC!clS`rP?9oY+%E4%iL0pdTleGjPaN{>;$x!7LFs2JO_ATrBWwP5}UTThLBrP%dUmH z*FP4QfZd+x7wm?ls3hXcppr0Fl#$BRg1wpvKx(WsStX^bbSI7UH^=cRqD?12QqZXO zDRvT=*wVWx2&3(8mpTc;Qevnn1O<|s3J&t_S`xEemH?HuTnT2LnK1^kxThuPEl_l~ zaWFi^tG9gA=grbImi%>5fBp3hx~UFUnS^gJ1_saSRoxx>^j^$&ttPyUvu#}trIwL6 zgxWkk?;sfp7g?nkPFa|Jve2m5(wT-S?e*Z9W##O}Qy+{6mEgr^^dD(-1+Epky~NAH zxh!52R{jV>6CJ{505K&dt6>w;X={de;5!B7?H`Z0y`jug2NO^o3JvvzFQ?%Kh z6s7(mou{D=?zN9S=hX3u$$-!A@-}?~*+tL4ao3oXb2pP_e;aDZcG*g_zx$h^h$T|1 z*#qy4Rjxw6vk|^w;PkZQH~WimuPw57(z+LvO4t;5R=sb~`o#EHN6^^!j=R3K*K^9J zh=++EXnWi5E>a?|Pg6=nL}sMB2V&_*{>WG7y?5OXr5MJP=FyH+2opk-zb@`dRtM-J zjpHAz!s}$pTEpA!cK=0Cr^XLfC^szW7J3X$M$M#5*Q?T|Yp0T}y z;*6k4Zw?*bpt^77n86hJ+|IL;t%;3*2=Gfji7;B9nca<{tA`JtnsqDO%$!O{Pj@1D5>N!M$M;^>~Vvb?F3#%%&_pBr7D- zz9HrUBza9aU#RGw0S&An{QBqi-)zmw(5s$OBTFF{kDcYSk3OE zWA02vKk8Q;=^^HoN+ZpOgmuTVHxKba(}O6V#NFJ*N&~vlJoOcTM>&ztLlyjykpm&m|j(HjwLN=sW8tYrCE#Twi4@rpwP59FBxe0!FaRJi4iyYz-1*I3Ft z^qZ{&YyiQjifj7%r967q4TsL>WJnI?OT5pZHv*638&swGv_*($CT!Z$&vT7reyS04 z6C%lV#&A|k^4_c#GwYZ#xiprrxBVC^Wn?8TLauxgw)WwoOhb0p{QnGg;D$%P6m(D{M~fRY~npNbC+}N zsNwc0?bOZWat?P87nyJG74B-lA2Rw5+fP3Z?h>jCtGo`{FL~DqVo7p^&Mi!;&^YyZ z9^y>*$~9AH8orP69+dJhb%}~j3Qf5^%)Yge>^K%_tS>9-4X~u?IoHhUJz@$E5K9D# zy7uI}XGZbXA}5TJ?l9hnvD4kjD+ zP*(*4abA|g5#mQ#xCZ229`U9Pxd&wEZA)oh={0E1ousuu71I=pe30*`OZk^j z>U-XWSfRtmJ?F&$XD}Y^572AXr)sgZ5~EKgVNee!3ExKzFBe<1SNao60Qpy%ervMf zSz?jIvq0}q4RHL$n+*-p#Od1MU6QB6v?=gLNdZt7PiT;sUm!E61GYujQVe6la%HKV zKY5p0$+@O)*vPx4zmc=OEOA0w^++5;Yb01|t-hE^Khy8#Z?m8}MhDo7y49KQJ^-9f z#@lRj>s0&%xzNfIa{c%d)^iffu;JxWYBgCN7RQd+4xE->S*jiAlK*7LJmM3txjv%J zMYYyn#7@Scjw25FfkQHn!#-!=(-ebdcB)DcvXc=Ys-Swc79+ApLv;-0GvuJomaPw^ zSZn7O#GC#4?Z4||rY>QOS~%~t9Zk4hSZ!z*R-A$9vCom;+)fP$sts?tZxD&okZf>m z#8BGrZ>9S;qf=7NIs|!APINM01MKsijHoD%F)SS0z!Yac8twllOU^)jIsCg|@R|Io+t5?by zx#6fIwuPE9rFFpIv;F|l>fL)A&PRm@%Ev>bkhUPl6(-dJi^gWqfP-QRR&0C3oaU$L z73D5GL6xN(O&q#jtl2}5_wE@U1XtCA*wX~CyB)O%K{`;q1bnUGXr(K;@Qq~dn|-W; zqOO}4hlbiH!ZOW4BBACM>Q)W!fRJ)rt<=~?4x+MZPlBb^mBWLH%U1OstJf&@`o;e_ z9IJRXsQ+pet%w8qqJA%fQIwP&ikeyemOVC7`bV?mAUv0ap52JFr`apZ2Ch-9qKOb%$ z3iG=Ne!h}$kgy`I_gR=WYOTnhB19Bs^O#Ath8KX*AmB3XM{9sL&WoOBq094q1N zfoI{rLIgA&H$UVFvOGN>xrxxPzIL69rWs8qHNO$Pyztc{2@+NzY7NG zE=QGzg^U^WM>suf!i0y+5~;Sq(HvxhpL7agmJWu%s=Yxj8sZ`km=X&&UDGT!HZp8O z^)LIW4dC9>d=K;V8(T*WKzygGmS$4Ja=4fp!h8I~azSd<2HJ_D2qZxc%uE?>8}uWD zWd{R~fzi4l98KA2*}Tld2tOp1fZMQfoe5bfKN)I{H$;u0v`)HTmqa8F>grO$Q}sOq zqYcljv*U1|R(#UV`BsVs&Y0C@HQWMhA6kI8--*GsGm=}f1$dg7(9fB8;k&sSb5K&l zL7-}~u~9Fo40GIyG~f#CHghkVwP-wKDB)6#`Q-2lm6^^4)^s`?v2|tx!8|y&p)5I4 z7qkQ=0F`Rkg>1WvsR?7VQ+n;qNQRoxkksg(@6q_Zq*ot0WZTaq2*zCICR`K7s(6h{ zHQ=hToFRGW^7_b{-q#y&(=}$l8UaD~Y<97eW`g2@tj33&3`d2GtILQr4eq}0-6Vb* z;IQ>o5y?KmxW%0E)9@7Bf7?p9PJw4fU-g76HX)VL7N9v#`DJ&^|tMk7#~1>f4R9| zqe;gHx6D~w%M@O8j^@|FC4?ee?BXd=T`ip7zr)8fI}>axtFSQRn}rY2DcDtAwOUBj?8m-?Kv% zYQOOU!C@c-LAH#U`m8B8;%)0e&3;Y)Z6>jr?Y4ff=%!>}N7Dl-y~BsNh8VmlMf&PO z6i@**xIN4)Dkx?_KHG)S$BenN@O2&iZ*gZ2@zhiJEJ&({!LdAQxWwX*!R4xTtD77YrkI4>a6n~&t?Bh&8oo|v$pF{QYVc#W`sLz<&1JD1rugN6&gbkc;BpHos`qzL#SekPw~-V7!Rpo znRT*=U^8MR=24DqCLgewUQ5bpoY|LhynQ4eqjIOPO{E+g@D8%IDa7Wc8f~Bn-tbb^4teCAy3sT zzZ%hRILEI)Gaa7>sdYa?X0doYI zc)6M*C!7!QQ3x9p1PzBxk3G#CIR z80;)?7$e72e{uF`;(ux1(;{A&CdSpqu2Jy@oOp36VQ*eo3<_Ha>$s* zHqM22G6nVeB)K|xk2;rN7eB7Y3X@*z4k0IX10Ej;<&PsRl}y`BA0 z?>dTKWk+>9FSh0vHiRxC`&{?1>DR>qHI z@cbo9?Os&*yh1IgVi6#&nVsATcf2fP{v0!xpG8a{N(@gn$U0t$<=RBkgQ)>ROGv4Q z@LdLxSTWmJO!SQjA-DZ;$n5*8SaA3v|7s14$B#EmT5vI#g|7>vdRSd7o4DQ!h9gdP zwU_a#Tm;k(pFzF7e~a7tAg{orI>tL~;o)f`t9{+r`*nYz*Hd7wJ-7b-wT^}&5qP<( z;xYBqfbbSnbhD5pHLs7=lx;~rUKhdCI4M}=O1j8l+lj>V8hZsBmi8cvuE&V!hy zZW`N~BUEFF$OCtsRdeh^$y|3Jv_AfwbBxHGawSVKw&^s$qz6n3ZALvAX%7U<3%)k# zDb}@|3Z^ENM_4W~4^TE5EP4vZltMR;DmcPDcNxgr$aF)Kmz1Wvjr9E}eY4`HT|15#asrQ`zgIw9Y z{N%m=D1plv_pba0`DOR=WB2}}1TIg!_sf5f?%m5{@BK#!T+X`J{r?j8lOi$GreoI4 z92pqSQzYhGH$$qfyuio6lQK$ah7QzzRet^Q>Q6yiH~mx1y2^Vmu6`f%#ioC%c~E(;=jx9^ z>o@&V&B@Ap-B))9E#LG{H4~NhUS16j^4;{m&U}2&`Eu}gz8_Bf&kf-NbLXsix# z3&Z_Q;&sIyQZystP2vs3UQ$qf=Dmg={-f;ZS#UM(^4ZjTu3sxxOnkZr_e1@Eq3dVw z-g8}}{AObJ9^8-hzw~5Zjr+eS|7+#f6KnV2_SEm`*?Tqa>e*}mjgIfZ1=Szw>AD(s z^K8Yv{58sr6aR&F^=!Htcjav6y|Ay9t0tE0srjM)ztD}d#rMM2DA!LM*i-Xk{r;ZP zt8xDq4Sua$JF$IF&7S(bJ+W8guAR;QZtR+lv3W=ttH5 zFNIg>H(HCoCauw~&j{Tk`Uj(pq+{CknW2Y79jg6HgpcVr+KT_@X!o%#>ob2jl+|(c z*Clr!FTK(B_s0K4Xjj(yqrW=cy|i?H>))$(FW&Mn{g?kavi!PBv|hE}NqC9AzqNQ3 zX|eWUdg%Wq^pL*4t$00YzxLs!(EXx^s{Ms~cwclBIi~hQ`cr9t0_TfJl z)ze?M6?>8rwGT5xBSi01`(1?pRP^kYVMWY9YRO5|first child +ep1<-with(demogr,epdata(start=dbirth, event=dch1, rcensor=dint, + birth=dbirth,id=id, + addvar=subset(demogr,select=c(-id,-dbirth)))) + +# creates a new episode-data structure with a time-varying factor +# variable relating to the status "never married" (not_marr) or +# "ever married" (marr) +ep2<-splitter(ep1,split=ep1$d1marr,tvar.lev=c("not_marr","marr"), + tvar.name="mar") + +# Estimates age profiles for the transition to the first birth +# according to the following factors: +# sex (respondent'sex w/2 levels: 'Male', 'Female'); +# edu ('Level of education w/3 levels: 'low_sec','upp_sec', 'tert'); +# mar (ever married w/2 levels: 'not_marr', 'marr') + +ch1.ap<-ageprofile(formula=~sex+edu+mar, epdata=ep2, + tr.name="First child", agelimits=c(15,50)) + +# Plot age profiles in three different graphs +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("Male","Female"),title='first child by sex') +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("low_sec","upp_sec","tert"),title='first child by education') +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("not_marr","marr"),title='first child by marital status', + ylim=0.4) +} +\keyword{package} +\keyword{smooth} +\keyword{age profiles} +\keyword{transition rates} +\keyword{GAM models} diff --git a/man/ageprofile.Rd b/man/ageprofile.Rd new file mode 100755 index 0000000..03f5cd5 --- /dev/null +++ b/man/ageprofile.Rd @@ -0,0 +1,174 @@ +\name{ageprofile} +\alias{ageprofile} +\title{Computes smoothed transition rates by respondent's age (age profiles).} +\description{ +Computes age profiles for a specific transition between two states according to +a set of time-fixed or time-varying covariates. It needs an micro-level episode-data +structure, i.e. a longitudinal dataset containing the following variables:\cr \cr +start : starting date of observation;\cr +Tstop : ending date of observation;\cr +Status : is 0 if right censored; 1 if event occurred; 2 if left censored; + 3 if interval censored. Status is equal to 1 if and only if the date of the event + precedes the date of right censoring.\cr +Agestart - ages at starting date (Start)\cr +Agestop Age at the ending date (Stop)\cr +An episode-data structure can be obtained through the command epdata.} +\usage{ +ageprofile(formula, epdata, tr.name = "Transition", win, + method = "car", agelimits = c(0, 100), + outfile = FALSE, tails = c(FALSE, FALSE), + subset = TRUE, weight, + sig.eff = TRUE, sig.lev = 0.05) +} +\arguments{ + \item{formula}{ +a formula object specifying the ~ operator and factors variables +separated by '+' operators. Response variable on the left is not required. +The expression ~1 implies that only the baseline (age profile for the whole sample) +will be provided. Note that interaction terms are not recognized by the formula: +~v1:v2 or ~v1*v2 gives the same results as ~v1+v2. +} + \item{epdata}{ +a dataframe containing episode-data prepared by epdata command optionally with +time-varying factor variables created by splitter command. +} + \item{tr.name}{a string containing the name of the considered transition} + \item{win}{ +a matrix with two columns containing the initial and final calendar date specifiyng +a restricted window of observation. Only events and exposure times referring to +this window will be considered in the analysis. If win argument is not specified, +the whole episode is considered. +For example, let us consider an episode data structure coming from a retrospective +survey held on January 1, 2010. Since we want to compute transition rates according to +behaviours experienced in the last 10 years, we can restrict our window ob observation +to the decade January 1, 2000 - January 1, 2010. Thus, the win argument would be a +matrix with two columns, the first containing the exact data 2000 and the second the exact +date 2010.\cr +Window of observation may be also limited to specific ages or events. For example, +win=cbind(date_at_birth+15,date_at_birth+20) restricts observations within the +age group 15-19 (completed) years of age and\cr +win=cbind(date_at_marriage,date_at_the_interview) restricts observatins +between marriage and the interview. + +Note that in the win matrix no missing data are allowed. +} + \item{method}{specifies the type of rates to compute. There are three + alternatives:\cr +'cpr' for Cohort Period Rates +'car' for Cohort Age Rates +'cap' for Cohort Age Probabilities +} + \item{agelimits}{ + a couple of values indicating the lower and the upper limit of age interval + to be considered. It can be useful to restrict age interval when the transition + is never or almost never experienced outside a specific age interval. For example, + the transition rates for the first child birth may be limited to the age interval + (15-50).} + \item{outfile}{ + if TRUE writes the output on a file named 'trname.txt' where 'trname' is the + string specified in the tr.name argument. +} + \item{tails}{ + a vector of two logical elements indicating respectively if the left and the + right tails must be flattened. This option may be useful if we can assume + that the transition rate is approximately zero at the borders of the considered + age interval. +} + \item{subset}{ + an optional vector specifying a subset of observations to be used in the + estimation process. +} + \item{weight}{ + an optional vector containing (post-sampling) weights. +} + \item{sig.eff}{ + if TRUE the age profile for a specific subgroup defined by factor levels is + fixed as identical to the baseline if the relative risk (for aech specific + age-subinterval) is not statistically significant at the level specified by + sig.lev. In other words, the relative risk in a specific age subinterval is zero + if the pvalue is higher than sig.lev (see details). +} + \item{sig.lev}{ + specifies the maximum level significance under which the relative risk in a + specific age-subinterval is non-zero (if sig.eff=TRUE). The default value is 0.05. +} +} +\details{ +p-values for the null hypothesis that the corresponding parameter is zero is +calculated with reference to the t distribution with the estimated residual +degrees of freedom for the model fit if the dispersion parameter has been estimated, +and the standard normal if not. +} +\value{ +Gives a list of objects: +\item{profiles}{a matrix containing the smoothed age profiles for the whole sample (baseline) +and for each factor level considered} +\item{unsmoothed}{a matrix containing the unsmoothed transition rates, i.e. the ratio +between occurences and exposures for each age and factor level} +\item{knot}{a matrix with the two knots (in column) for each factor level (rows)} +\item{event}{a matrix with the number of event occured in the three age sub-intervals defined +by knots (columns) for each factor level (rows)} +\item{rrisk}{a matrix with the estimated relative risks in the three age sub-intervals + (columns) for each factor level (rows)} +\item{se}{a matrix with the standard error related to the relative risk computed +in the three age sub-intervals (columns) for each factor level (rows)} +\item{pvalue}{a matrix with the pvalues related to the hypothesis that the +relative risk in each age sub-interval is different from zero (columns) for each +factor level (rows)} +\item{factors}{a vector of names containing the factors specified in the formula argument} +\item{trname}{the name of the transition in analysisi as specified in the tr.name argument.} +\item{ANOVAtest}{pvalues anova test for each factor covariate} +\item{factors}{vector of string containing the names of factor covariates as specified in +the formula argument.} +\item{method}{the method used for the computaion of rates as speficied in the method argument.} +\item{agelimits}{the age interval as specified in the agelimits argument.} +\item{tails}{tails argument.} +} +\references{Impicciatore R. and Billari F.C., (2010), MAPLES: a general method for the estimation +of age profiles from standard demographic surveys (with an application to fertility),DEAS WP, \cr +http://ideas.repec.org/p/mil/wpdepa/2010-40.html} +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} +\note{ +In order to have an unambiguosly output, it is strongly raccomended to label +each level of factor variables and to avoid using the same label for different +factors. +} + +\seealso{ +\code{\link{epdata}}, +\code{\link{splitter}}, +\code{\link{plotap}} +} +\examples{ +# creates an episode-data structure relating to the transition +# childless-->first child +ep1<-with(demogr,epdata(start=dbirth, event=dch1, rcensor=dint, + birth=dbirth,id=id, + addvar=subset(demogr,select=c(-id,-dbirth)))) + +# creates a new episode-data structure with a time-varying factor +# variable relating to the status "never married" (not_marr) or +# "ever married"(marr) +ep2<-splitter(ep1,split=ep1$d1marr,tvar.lev=c("not_marr","marr"), + tvar.name="mar") + +# Estimates age profiles for the transition to the first birth +# according to the following factors: +# sex (respondent'sex w/2 levels: 'Male', 'Female'); +# edu ('Level of education w/3 levels: 'low_sec','upp_sec', 'tert'); +# mar (ever married w/2 levels: 'not_marr', 'marr') +ch1.ap<-ageprofile(formula=~sex+edu+mar, epdata=ep2, + tr.name="First child", agelimits=c(15,50)) + +# The estimates are obtained under the hypothesis of independence among +# factors. We can relax this hp by considering the interaction between +# factors. The following commands add the interaction between sex and edu. +ep2$inter<-ep2$sex:ep2$edu +ch1.ap<-ageprofile(formula=~sex+edu+mar+inter, epdata=ep2, + tr.name="First child", agelimits=c(15,50)) +} +\keyword{smooth} +\keyword{age|profile} +\keyword{transition|rates} diff --git a/man/demogr.Rd b/man/demogr.Rd new file mode 100755 index 0000000..3225243 --- /dev/null +++ b/man/demogr.Rd @@ -0,0 +1,59 @@ +\name{demogr.RData} +\alias{demogr.RData} +\alias{demogr} +\docType{data} +\title{Longitudinal data on marriage and first child birth.} +\description{ +Longitudinal data on marriage and first child birth. Individuals interwieved +in March 2004. +} +\usage{demogr} +\format{ +A data frame with 2017 observations on 8 variables. + +[,1] id integer ID \cr +[,2] weight numeric Individual weight \cr +[,3] dbirth numeric respondent's date at birth \cr +[,4] dint numeric date at the interview \cr +[,5] d1marr numeric date at marriage \cr +[,6] dch1 numeric date at first child birth \cr +[,7] sex factor sex \cr +[,8] edu factor level od efucation \cr + +} +\details{ +Individuals were interviewed at the end of 2003 retrospectively on family and +fertility life trajectory. \cr + +id: individual identification number (ID) \cr +weight: individual post-sampling weight (with mean = 1) \cr +dbirth: respondent's date at birth in exact years since January 1, 1900 \cr +dint: date at the interview in exact years since January 1, 1900 (March 2004 +for all respondents) \cr +d1marr: date at marriage (if any) in exact years since January 1, 1900. If NA, +the individual has never been married before the interview. \cr +dch1: date at marriage (if any) in exact years since January 1, 1900. If NA, +the individual is childless at the interview. \cr +sex: respondent's sex (factor w/2 levels: "Male", "Female") \cr +edu: respondent's level of education (factor w/3 levels: "low_sec","upp_sec","tert") \cr + +Dates in exact years have been computed considering the midpoint of a specific month. +Thus, March 1995 means March 15, 1995 and the date in exact years is 95+2.5/12=95.21. +} +\source{ +Synthetic data based on the format of the Generation and Gender Survey (GGS) (Vikat et al. 2007). +} +\references{ +Vikat A., Speder Z., Beets G., Billari F.C., Buhler C., Desesquelles A., Fokkema T., +Hoem J., MacDonald A., Neyer G., Pailhe A., Pinnelli A., Solaz A. (2007). +"Generations and Gender Survey (GGS): Toward a better understanding of relationship +and processes in the life course", Demographic research, 17 (14): 389-440. +\cr +\cr +Impicciatore R. and Billari F.C., (2010), MAPLES: a general method for the estimation +of age profiles from standard demographic surveys (with an application to fertility), DEAS WP,\cr +http://ideas.repec.org/p/mil/wpdepa/2010-40.html} +\examples{ +str(demogr) +} +\keyword{datasets} diff --git a/man/epdata.Rd b/man/epdata.Rd new file mode 100755 index 0000000..6c6789a --- /dev/null +++ b/man/epdata.Rd @@ -0,0 +1,72 @@ +\name{epdata} +\alias{epdata} +\title{Prepares episode data for event history analysis.} +\description{Starting from a micro-level dataset containing dates for some +relevant events, it prepares a file with an episode-data format suitable for +event history analysis.} +\usage{epdata(start, event, lcensor, rcensor, + subset = TRUE, birth, id, addvar) +} +\arguments{ + \item{start}{a vector of dates at which the unit starts to be at risk to experience + the transition. It may be unknown in case of left-censoring. When the episode + is no left censored, start coincides with the beginning of the episode.} + \item{event}{a vector of containing the date of the event of interest. It may + be unknown in case of right-censoring.} + \item{lcensor}{an optional vector of left censoring dates.} + \item{rcensor}{an optional vector of right censoring dates.} + \item{subset}{a logical expression indicating elements or rows to keep. + It specifies the subset of cases at risk of experiencing the event of interest. + If TRUE, the whole sample will be considered.} + \item{birth}{a vector containing the respondent date of birth. + This information in useful in order to compute ages at various events (NA not allowed).} + \item{id}{a vector of identification numbers. If id remains unspecified, a + default id is implemented.} + \item{addvar}{a dataframe containing additional variables to add in the resulting output.} +} +\details{All the specified objects must have the same length (i.e. the number of + cases must be the same). For any observation at least one between start and + lcensor and one between event and rcensor must be known. In general, we have + start<=lcensor<= event<=rcensor.} +\value{A data.frame suitable for MAPLES and R survival package. Following the + 'counting process' formulation of survival analysis, each record ha the + following variables: +\item{Tstart}{episode starting date} +\item{Tstop}{episode ending date} +\item{status}{0: right censored: 1: event occurred: 2: left censored: 3: interval censored. + Is equal to 1 if and only if the date of the event precedes the date of right censoring.} +\item{Agestart}{Age at starting date (Tstart)} +\item{Agestop}{Age at the ending date (Tstop)} +} +\references{Impicciatore R. and Billari F.C., (2010), MAPLES: a general method for the +estimation of age profiles from standard demographic surveys (with an application to fertility), +DEAS WP, \cr +http://ideas.repec.org/p/mil/wpdepa/2010-40.html} +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} +\seealso{ +\code{\link{splitter}}, +\code{\link{ageprofile}}, +\code{\link{plotap}} +} +\examples{ +# creates an episode-data structure relating +# to the transition childless-->first child +ep1<-epdata(start=demogr$dbirth, + event=demogr$dch1, + rcensor=demogr$dint, + birth=demogr$dbirth, + id=demogr$id, + addvar=subset(demogr,select=c(-id,-dbirth))) + +# or, using 'with' +ep1<-with(demogr,epdata(start=dbirth, event=dch1, + rcensor=dint, birth=dbirth,id=id, + addvar=subset(demogr,select=c(-id,-dbirth)))) +} +\keyword{datagen} +\keyword{episode data} +\keyword{event history analysis} +\keyword{survival} + diff --git a/man/listvar.Rd b/man/listvar.Rd new file mode 100755 index 0000000..3522d0c --- /dev/null +++ b/man/listvar.Rd @@ -0,0 +1,21 @@ +\name{listvar} +\alias{listvar} +\title{Shows variables in a dataframe.} +\description{Shows variables in a dataframe and the related number of column.} +\usage{listvar(df)} + +\arguments{ + \item{df}{is the dataframe.} +} + +\value{ +a dataframe with two columns: VAR containing the names of the variables and COL +with the number of related column in the original dataframe. +} + +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} +\examples{ +listvar(demogr) +} diff --git a/man/mkdate.Rd b/man/mkdate.Rd new file mode 100755 index 0000000..8ff7f1a --- /dev/null +++ b/man/mkdate.Rd @@ -0,0 +1,36 @@ +\name{mkdate} +\alias{mkdate} +\title{Compute dates.} +\description{ +Computes dates in continuous years or CMC (Century Month Code) since January 1, 1900. +} +\usage{mkdate(year, month, cmc = FALSE)} +\arguments{ + \item{year}{ +a vector containing the year of a date(XXXX format) +} + \item{month}{ +a vector containing the month of a date (XX format) +} + \item{cmc}{ +if TRUE date in CMC from January 1, 1900 is computed. Otherwise, date in continuos time +since January 1, 1900 is computed. +} +} + +\value{a vector containing dates in continuous time or CMC.} +\references{Impicciatore R. and Billari F.C., (2010), MAPLES: a general method for the estimation +of age profiles from standard demographic surveys (with an application to fertility), +DEAS WP, \cr +http://ideas.repec.org/p/mil/wpdepa/2010-40.html} +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} + +\examples{ +dy<-c(1996,2001,2005,2003) +dm<-c(3,9,5,12) +mkdate(year=dy,month=dm) +mkdate(year=dy,month=dm, cmc=TRUE) +} +\keyword{dates} diff --git a/man/plotap.Rd b/man/plotap.Rd new file mode 100755 index 0000000..df9078a --- /dev/null +++ b/man/plotap.Rd @@ -0,0 +1,109 @@ +\name{plotap} +\alias{plotap} +\title{Plots age profiles.} +\description{ +Plots age profiles estimated through command ageprofile. +} +\usage{ +plotap(x, lev.labels, baseline = TRUE, unsmoothed = FALSE, + xlim, ylim, title) +} + +\arguments{ + \item{x}{the resulting list given by the ageprofile command.} + \item{lev.labels}{ +a vector of strings specifying the factor levels to consider. For example: \cr +lev.labels=c("Male","Female","low_sec","upp_sec","tert") \cr +specifies that 5 curves will bedrawn, one for each level specified. +It is also possible to draw a combination of levels (under the hypothesis of +independence between factors)by inserting the symbol "*" between two or more levels. +For example: the string "Male*low_sec" draws a curve for the subgroups of men with +a lower secondary level of education; the string "Male*low_sec*worker" draws a curve for +the subgroup of lower educated men currently employed as a manual worker. +Even though it is possible to draw any number of curves in one graph, we raccomend +to consider no more than 5 levels (or combination of levels) in one graph.} + \item{baseline}{if TRUE the baseline will be drawn.} + \item{unsmoothed}{if TRUE the unsmoothed transition rates, i.e. the ratio +between occurences and exposures for each age and factor level, are plotted as points +in the graph.} + \item{xlim}{ +a vector of two values defining the limits of X axis measured in years of age +(default value: c(min(age), max(age))} + \item{ylim}{ +a value defining the upper limit of the Y axis (transition rates). +Given that transition rates a re strictly positive, the lower limit is always 0 +(default value: 2*max(baseline rates)). } + \item{title}{ +title of the graph (default: transition name as stored in the ageprofile argument. } +} +\references{Impicciatore R. and Billari F.C., (2010), MAPLES: a general method for the estimation +of age profiles from standard demographic surveys (with an application to fertility), +DEAS WP, \cr +http://ideas.repec.org/p/mil/wpdepa/2010-40.html} +\author{ +Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} +\seealso{ +\code{\link{epdata}}, +\code{\link{splitter}}, +\code{\link{ageprofile}}, +} +\examples{ +# creates an episode-data structure relating to the +# transition childless-->first child +ep1<-with(demogr,epdata(start=dbirth, event=dch1, rcensor=dint, + birth=dbirth,id=id, + addvar=subset(demogr,select=c(-id,-dbirth)))) + +# creates a new episode-data structure with a time-varying factor variable +# relating to the status "never married"(not_marr) or "ever married"(marr) +ep2<-splitter(ep1,split=ep1$d1marr,tvar.lev=c("not_marr","marr"), + tvar.name="mar") + +# Estimates age profiles for the transition to the first birth +# according to the following factors: +# sex (respondent'sex w/2 levels: 'Male', 'Female'); +# edu ('Level of education w/3 levels: 'low_sec','upp_sec', 'tert'); +# mar (ever married w/2 levels: 'not_marr', 'marr') + +ch1.ap<-ageprofile(formula=~sex+edu+mar, epdata=ep2, + tr.name="First child", agelimits=c(15,50)) + +# Plot age profiles in three different graphs +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("Male","Female"),title='first child by sex') +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("low_sec","upp_sec","tert"),title='first child by education') +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("not_marr","marr"),title='first child by marital status', + ylim=0.4) + +# Plot age profiles for the combined effect of sex and level of education +# under the independence hypothesis +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("Female*low_sec","Female*upp_sec","Female*tert"), + title='first child by education - women (indep hp)') +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("Male*low_sec","Male*upp_sec","Male*tert"), + title='first child by education - men (indep hp)') + +# The estimates are obtained under the hypothesis of independence among +# factors. We can relax this hp by considering the interaction between +# factors. The following commands add the interaction between sex and edu. +ep2$inter<-ep2$sex:ep2$edu +ch1.ap<-ageprofile(formula=~sex+edu+mar+inter, epdata=ep2, + tr.name="First child", agelimits=c(15,50)) + +# Plot age profile for the interaction between sex and level of education +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("Female:low_sec","Female:upp_sec","Female:tert"), + title='first child by education - women') +plotap(ch1.ap,base=TRUE, unsmoo=TRUE, + lev=c("Male:low_sec","Male:upp_sec","Male:tert"), + title='first child by education - men') + +} +\keyword{smooth} +\keyword{age|profiles} +\keyword{transition|rates} diff --git a/man/splitter.Rd b/man/splitter.Rd new file mode 100755 index 0000000..d0f7a81 --- /dev/null +++ b/man/splitter.Rd @@ -0,0 +1,75 @@ +\name{splitter} +\alias{splitter} +\title{Creates a time-varying factor variables within a episode-data.} +\description{ +Splits episode in sub-episodes according to the dates specified as input. + The resulting dataframe contains a new time-varying factor variable with different + user-defined values in each subinterval. +} +\usage{splitter(epdata, split, tvar.name = "Tvar", tvar.lev)} +\arguments{ + \item{epdata}{a dataframe containing an episode-data file (obtained through the epdata command)} + \item{split}{an array containing the date(s) at one (or more) specific event(s) + that causes the change in the resulting time-varying variable. The original episode + will be split at this date(s). If n vectors of dates have been spcified, the original + episode will be splitted in n+1 subepisodes. For example, let us consider a matrix + with three columns containing the dates at the first, second, and third child birth. + The original episode may be split into (up to) four subepisodes according to non-NA + specified dates (see details).} + \item{tvar.name}{a string containing the name of the time-varying factor + variable in the resulting dataframe.} + \item{tvar.lev}{a vector of n+1 strings containing the level labels of the resulting + time-varying variable (where n is the number of columns in the split argument).} +} +\details{ +Split dates refer to event that can be experienced sequentially (first birth, +second birth, third birth, ecc or first marriage, second marriage, third marriage, etc.). +This implies that \cr +1. Split dates must be strictly sequential: date1<=date2<=date3<=etc. \cr +2. NA date cannot be followed by a non-NA date (e.g. for the j-th the sequence +of dates "date1,date2, NA" is allowed and it means that event3 has not been +experienced, whereas "NA, date1, date2" is not allowed because we do not have +any information about the first change ). \cr +The number of rows in the split argument must be the same as the number of rows +in the epdata argument. +} +\value{ +In the resulting dataframe each row is a subepisode obtained through the splitting +procedure. The columns will be the same as the epdata argument plus one time-varying +factor variable with levels specified through tvar.lev argument. \cr +In order to create several time-varying variables on the same dataframe, splitter +procedure should be applied repeatedly for each new time-varying variable. At any +step the input dataframe will be the resulting dataframe in the previous one. +} +\references{Impicciatore R. and Billari F.C., (2010), MAPLES: a general method for the estimation +of age profiles from standard demographic surveys (with an application to fertility), +DEAS WP, \cr +http://ideas.repec.org/p/mil/wpdepa/2010-40.html} + +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} +\seealso{ +\code{\link{epdata}}, +\code{\link{ageprofile}}, +\code{\link{plotap}} +} + +\examples{ +# creates an episode-data structure relating to the transition +# childless-->first child +ep1<-with(demogr,epdata(start=dbirth, event=dch1, rcensor=dint, + birth=dbirth,id=id, + addvar=subset(demogr,select=c(-id,-dbirth)))) + +# creates a new episode-data structure with a time-varying factor +# variable relating to the status "never married" (not_marr) or +# "ever married" (marr) +ep2<-splitter(ep1,split=ep1$d1marr,tvar.lev=c("not_marr","marr"), + tvar.name="mar") +} +\keyword{datagen} +\keyword{episode|data} +\keyword{event|history|analysis} +\keyword{survival} +\keyword{episode|splitting} diff --git a/man/tabm.Rd b/man/tabm.Rd new file mode 100755 index 0000000..b7bcb32 --- /dev/null +++ b/man/tabm.Rd @@ -0,0 +1,31 @@ +\name{tabm} +\alias{tabm} +\title{Print regression estimates for previously fitted linear and logit +regression models.} +\description{ +For linear regression models, it prints some basic information about the model +estimates ( Beta, pvalues, AIC, adjusted R square); for logit models, it also prints +exp(Beta).} +\usage{tabm(mod, pvalues = TRUE, digits = 3)} +\arguments{ + \item{mod}{is an oject of class 'lm' or 'glm'.} + \item{pvalues}{if TRUE pvalues will be printed out.} + \item{digits}{specifies the number of digits in the output.} +} +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} + +\seealso{ +\code{\link{summary.lm}}, +\code{\link{summary.glm}}, +} + +\examples{ +chi<-ifelse(!is.na(demogr$dch1),1,0) +logit.mod<- glm(chi ~ demogr$sex + demogr$edu, + family=binomial(link = "logit")) +tabm(logit.mod) +# for more details +summary(logit.mod) +} diff --git a/man/tabx.Rd b/man/tabx.Rd new file mode 100755 index 0000000..3f4360a --- /dev/null +++ b/man/tabx.Rd @@ -0,0 +1,28 @@ +\name{tabx} +\alias{tabx} +\title{Prints univariate and bivariate frequency table.} +\description{Prints univariate or a bivariate frequency distribution table including marginal +distribution and total number of cases.} +\usage{tabx(x, y, prow = FALSE, pcol = FALSE, chisq=FALSE )} +\arguments{ + \item{x}{a vector which can be interpred as factor.} + \item{y}{an optional second vector which can be interpreted as factor (columns + in the contingency table).} + \item{prow}{if TRUE it adds a bivariate table containing row percentages.} + \item{pcol}{if TRUE it adds a bivariate table containing column percentages.} + \item{chisq}{gives the results for a Chi square test for independence of all factors.} +} +\author{Roberto Impicciatore +\email{roberto.impicciatore@unimi.it} +} + +\seealso{ +\code{\link{table}}, +\code{\link{summary.table}}, +} +\examples{ +tabx(demogr$sex,demogr$edu,chisq=TRUE) +tabx(demogr$sex,demogr$edu,prow=TRUE,pcol=TRUE) +} +\keyword{univar} +\keyword{frquencies}