R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("gss-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('gss') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "LakeAcid" > > ### * LakeAcid > > flush(stderr()); flush(stdout()) > > ### Name: LakeAcidity > ### Title: Water Acidity in Lakes > ### Aliases: LakeAcidity > ### Keywords: datasets > > ### ** Examples > > ## Converting latitude and longitude to x-y coordinates > ## Not run: > ##D convert <- function(lat, lon) { > ##D lat <- lat/180*pi > ##D lon <- lon/180*pi > ##D m.lat <- (max(lat)+min(lat))/2 > ##D m.lon <- (max(lon)+min(lon))/2 > ##D x <- cos(m.lat)*sin(m.lon-lon) > ##D y <- sin(lat-m.lat) > ##D cbind(x,y) > ##D } > ##D data(LakeAcidity) > ##D convert(LakeAcidity$lat,LakeAcidity$lon) > ##D ## Clean up > ##D rm(LakeAcidity,convert) > ## End(Not run) > > > > cleanEx(); ..nameEx <- "gssanova" > > ### * gssanova > > flush(stderr()); flush(stdout()) > > ### Name: gssanova > ### Title: Fitting Smoothing Spline ANOVA Models with Non Gaussian > ### Responses > ### Aliases: gssanova > ### Keywords: models regression smooth > > ### ** Examples > > ## Fit a cubic smoothing spline logistic regression model > test <- function(x) + {.3*(1e6*(x^11*(1-x)^6)+1e4*(x^3*(1-x)^10))-2} > x <- (0:100)/100 > p <- 1-1/(1+exp(test(x))) > y <- rbinom(x,3,p) > logit.fit <- gssanova(cbind(y,3-y)~x,family="binomial") > ## The same fit > logit.fit1 <- gssanova(y/3~x,"binomial",weights=rep(3,101)) > ## Obtain estimates and standard errors on a grid > est <- predict(logit.fit,data.frame(x=x),se=TRUE) > ## Plot the fit and the Bayesian confidence intervals > plot(x,y/3,ylab="p") > lines(x,p,col=1) > lines(x,1-1/(1+exp(est$fit)),col=2) > lines(x,1-1/(1+exp(est$fit+1.96*est$se)),col=3) > lines(x,1-1/(1+exp(est$fit-1.96*est$se)),col=3) > ## Clean up > ## Not run: > ##D rm(test,x,p,y,logit.fit,logit.fit1,est) > ##D dev.off() > ## End(Not run) > > > > cleanEx(); ..nameEx <- "gssanova1" > > ### * gssanova1 > > flush(stderr()); flush(stdout()) > > ### Name: gssanova1 > ### Title: Fitting Smoothing Spline ANOVA Models with Non Gaussian > ### Responses > ### Aliases: gssanova1 > ### Keywords: models regression smooth > > ### ** Examples > > ## Fit a cubic smoothing spline logistic regression model > test <- function(x) + {.3*(1e6*(x^11*(1-x)^6)+1e4*(x^3*(1-x)^10))-2} > x <- (0:100)/100 > p <- 1-1/(1+exp(test(x))) > y <- rbinom(x,3,p) > logit.fit <- gssanova1(cbind(y,3-y)~x,family="binomial") > ## The same fit > logit.fit1 <- gssanova1(y/3~x,"binomial",weights=rep(3,101), + id.basis=logit.fit$id.basis) > ## Obtain estimates and standard errors on a grid > est <- predict(logit.fit,data.frame(x=x),se=TRUE) > ## Plot the fit and the Bayesian confidence intervals > plot(x,y/3,ylab="p") > lines(x,p,col=1) > lines(x,1-1/(1+exp(est$fit)),col=2) > lines(x,1-1/(1+exp(est$fit+1.96*est$se)),col=3) > lines(x,1-1/(1+exp(est$fit-1.96*est$se)),col=3) > > ## Fit a mixed-effect logistic model > data(bacteriuria) > bact.fit <- gssanova1(infect~trt+time,family="binomial",data=bacteriuria, + id.basis=(1:820)[bacteriuria$id%in%c(3,38)],random=~1|id) > ## Predict fixed effects > predict(bact.fit,data.frame(time=2:16,trt=as.factor(rep(1,15))),se=TRUE) $fit [1] -1.0412810 -0.8808566 -0.7274899 -0.5863526 -0.4626204 -0.3594481 [7] -0.2852688 -0.2422548 -0.2225406 -0.2170295 -0.2237757 -0.2380233 [13] -0.2527194 -0.2639604 -0.2750067 $se.fit [1] 0.2141301 0.1807036 0.1655871 0.1612462 0.1614710 0.1631558 0.1652667 [8] 0.1676805 0.1705044 0.1739185 0.1786438 0.1868461 0.2027299 0.2319145 [15] 0.2788945 > ## Estimated random effects > bact.fit$b [1] 0.053896069 -0.059707512 0.002650160 0.019476034 0.527924713 [6] 0.002650160 0.019476034 -0.236458817 -0.460964737 0.214477671 [11] 0.313803198 0.124697953 0.088919784 -0.018811134 -0.059707512 [16] -0.148392912 0.151390792 -0.624396817 0.002650160 0.133004907 [21] 0.375725062 -0.263136395 0.002650160 0.002650160 0.053896069 [26] -0.148392912 -0.067183018 0.133004907 0.088919784 -0.302622176 [31] -0.146956654 0.416798495 -0.086023593 0.141972821 -0.099048152 [36] -0.148832753 -0.048843831 0.809005657 0.407426955 0.063743872 [41] -0.052941540 0.443049542 -0.360126991 -0.078877298 0.621790900 [46] 0.064058918 0.408643909 0.621790900 -0.208249061 -0.480887807 [51] 0.106721060 -0.099620071 -0.256780447 -0.208249061 -0.052941540 [56] 0.443049542 -1.553417982 0.106721060 -0.480887807 -0.315591097 [61] -0.256780447 0.620441328 0.621790900 0.443049542 -0.020832185 [66] -0.950563051 0.063743872 -0.415538945 0.443049542 0.106721060 [71] -0.378771301 -0.174898092 > > ## Clean up > ## Not run: > ##D rm(test,x,p,y,logit.fit,logit.fit1,est,bacteriuria,bact.fit) > ##D dev.off() > ## End(Not run) > > > > cleanEx(); ..nameEx <- "mkran" > > ### * mkran > > flush(stderr()); flush(stdout()) > > ### Name: mkran > ### Title: Generating Random Effects in Mixed-Effect Models > ### Aliases: mkran > ### Keywords: models regression > > ### ** Examples > > ## Toy data > test <- data.frame(grp=as.factor(rep(1:2,c(2,3)))) > ## First formula > ran.test <- mkran(~1|grp,test) > ran.test$z [,1] [,2] [1,] 1 0 [2,] 1 0 [3,] 0 1 [4,] 0 1 [5,] 0 1 > ran.test$sigma$fun(2,ran.test$sigma$env) # diag(10^(-2),2) [,1] [,2] [1,] 0.01 0.00 [2,] 0.00 0.01 > ## Second formula > ran.test <- mkran(~grp|grp,test) > ran.test$z [,1] [,2] [1,] 1 0 [2,] 1 0 [3,] 0 1 [4,] 0 1 [5,] 0 1 > ran.test$sigma$fun(c(1,2),ran.test$sigma$env) # diag(10^(-1),10^(-2)) [,1] [,2] [1,] 0.1 0.00 [2,] 0.0 0.01 > ## Clean up > ## Not run: rm(test,ran.test) > > > > cleanEx(); ..nameEx <- "predict.ssanova" > > ### * predict.ssanova > > flush(stderr()); flush(stdout()) > > ### Name: predict.ssanova > ### Title: Predicting from Smoothing Spline ANOVA Fits > ### Aliases: predict.ssanova predict.ssanova1 > ### Keywords: models regression smooth > > ### ** Examples > > ## THE FOLLOWING EXAMPLE IS TIME-CONSUMING > ## Not run: > ##D ## Fit a model with thin-plate marginals, where geog is 2-D > ##D data(LakeAcidity) > ##D fit <- ssanova(ph~log(cal)*geog,"tp",LakeAcidity) > ##D ## Obtain estimates and standard errors on a grid > ##D new <- data.frame(cal=1,geog=I(matrix(0,1,2))) > ##D new <- model.frame(~log(cal)+geog,new) > ##D predict(fit,new,se=TRUE) > ##D ## Evaluate the geog main effect > ##D predict(fit,new,se=TRUE,inc="geog") > ##D ## Evaluate the sum of the geog main effect and the interaction > ##D predict(fit,new,se=TRUE,inc=c("geog","log(cal):geog")) > ##D ## Evaluate the geog main effect on a grid > ##D grid <- seq(-.04,.04,len=21) > ##D new <- model.frame(~geog,list(geog=cbind(rep(grid,21),rep(grid,rep(21,21))))) > ##D est <- predict(fit,new,se=TRUE,inc="geog") > ##D ## Plot the fit and standard error > ##D par(pty="s") > ##D contour(grid,grid,matrix(est$fit,21,21),col=1) > ##D contour(grid,grid,matrix(est$se,21,21),add=TRUE,col=2) > ##D ## Clean up > ##D rm(LakeAcidity,fit,new,grid,est) > ##D dev.off() > ## End(Not run) > > > > cleanEx(); ..nameEx <- "ssanova" > > ### * ssanova > > flush(stderr()); flush(stdout()) > > ### Name: ssanova > ### Title: Fitting Smoothing Spline ANOVA Models > ### Aliases: ssanova > ### Keywords: smooth models regression > > ### ** Examples > > ## Fit a cubic spline > x <- runif(100); y <- 5 + 3*sin(2*pi*x) + rnorm(x) > cubic.fit <- ssanova(y~x,method="m") > ## The same fit with different internal makes > tp.fit <- ssanova(y~x,"tp",method="m") > ## Obtain estimates and standard errors on a grid > new <- data.frame(x=seq(min(x),max(x),len=50)) > est <- predict(cubic.fit,new,se=TRUE) > ## Plot the fit and the Bayesian confidence intervals > plot(x,y,col=1); lines(new$x,est$fit,col=2) > lines(new$x,est$fit+1.96*est$se,col=3) > lines(new$x,est$fit-1.96*est$se,col=3) > ## Clean up > ## Not run: > ##D rm(x,y,cubic.fit,tp.fit,new,est) > ##D dev.off() > ## End(Not run) > > ## Fit a tensor product cubic spline > data(nox) > nox.fit <- ssanova(log10(nox)~comp*equi,data=nox) > ## Fit a spline with cubic and nominal marginals > nox$comp<-as.factor(nox$comp) > nox.fit.n <- ssanova(log10(nox)~comp*equi,data=nox) > ## Fit a spline with cubic and ordinal marginals > nox$comp<-as.ordered(nox$comp) > nox.fit.o <- ssanova(log10(nox)~comp*equi,data=nox) > ## Clean up > ## Not run: rm(nox,nox.fit,nox.fit.n,nox.fit.o) > > > > cleanEx(); ..nameEx <- "ssanova1" > > ### * ssanova1 > > flush(stderr()); flush(stdout()) > > ### Name: ssanova1 > ### Title: Fitting Smoothing Spline ANOVA Models > ### Aliases: ssanova1 > ### Keywords: smooth models regression > > ### ** Examples > > ## Fit a cubic spline > x <- runif(100); y <- 5 + 3*sin(2*pi*x) + rnorm(x) > cubic.fit <- ssanova1(y~x) > ## Obtain estimates and standard errors on a grid > new <- data.frame(x=seq(min(x),max(x),len=50)) > est <- predict(cubic.fit,new,se=TRUE) > ## Plot the fit and the Bayesian confidence intervals > plot(x,y,col=1); lines(new$x,est$fit,col=2) > lines(new$x,est$fit+1.96*est$se,col=3) > lines(new$x,est$fit-1.96*est$se,col=3) > ## Clean up > ## Not run: > ##D rm(x,y,cubic.fit,new,est) > ##D dev.off() > ## End(Not run) > > ## Fit a tensor product cubic spline > data(nox) > nox.fit <- ssanova1(log10(nox)~comp*equi,data=nox) > ## Fit a spline with cubic and nominal marginals > nox$comp<-as.factor(nox$comp) > nox.fit.n <- ssanova1(log10(nox)~comp*equi,data=nox) > ## Fit a spline with cubic and ordinal marginals > nox$comp<-as.ordered(nox$comp) > nox.fit.o <- ssanova1(log10(nox)~comp*equi,data=nox) > ## Clean up > ## Not run: rm(nox,nox.fit,nox.fit.n,nox.fit.o) > > > > cleanEx(); ..nameEx <- "ssden" > > ### * ssden > > flush(stderr()); flush(stdout()) > > ### Name: ssden > ### Title: Estimating Probability Density Using Smoothing Splines > ### Aliases: ssden > ### Keywords: smooth models distribution > > ### ** Examples > > ## 1-D estimate: Buffalo snowfall > data(buffalo) > buff.fit <- ssden(~buffalo,domain=data.frame(buffalo=c(0,150))) > plot(xx<-seq(0,150,len=101),dssden(buff.fit,xx),type="l") > plot(xx,pssden(buff.fit,xx),type="l") > plot(qq<-seq(0,1,len=51),qssden(buff.fit,qq),type="l") > ## Clean up > ## Not run: > ##D rm(buffalo,buff.fit,xx,qq) > ##D dev.off() > ## End(Not run) > > ## 2-D with triangular domain: AIDS incubation > data(aids) > ## rectangular quadrature > quad.pt <- expand.grid(incu=((1:40)-.5)/40*100,infe=((1:40)-.5)/40*100) > quad.pt <- quad.pt[quad.pt$incu<=quad.pt$infe,] > quad.wt <- rep(1,nrow(quad.pt)) > quad.wt[quad.pt$incu==quad.pt$infe] <- .5 > quad.wt <- quad.wt/sum(quad.wt)*5e3 > ## additive model (pre-truncation independence) > aids.fit <- ssden(~incu+infe,data=aids,subset=age>=60, + domain=data.frame(incu=c(0,100),infe=c(0,100)), + quad=list(pt=quad.pt,wt=quad.wt)) > ## conditional (marginal) density of infe > jk <- cdssden(aids.fit,xx<-seq(0,100,len=51),data.frame(incu=50)) > plot(xx,jk$pdf,type="l") > ## conditional (marginal) quantiles of infe (TIME-CONSUMING) > ## Not run: > ##D cqssden(aids.fit,c(.05,.25,.5,.75,.95),data.frame(incu=50),jk$int) > ## End(Not run) > ## Clean up > ## Not run: > ##D rm(aids,quad.pt,quad.wt,aids.fit,jk,xx) > ##D dev.off() > ## End(Not run) > > > > cleanEx(); ..nameEx <- "sshzd" > > ### * sshzd > > flush(stderr()); flush(stdout()) > > ### Name: sshzd > ### Title: Estimating Hazard Function Using Smoothing Splines > ### Aliases: sshzd > ### Keywords: smooth models survival > > ### ** Examples > > ## Model with interaction > data(gastric) > gastric.fit <- sshzd(Surv(futime,status)~futime*trt,data=gastric) > ## exp(-Lambda(600)), exp(-(Lambda(1200)-Lambda(600))), and exp(-Lambda(1200)) > survexp.sshzd(gastric.fit,c(600,1200,1200),data.frame(trt=as.factor(1)),c(0,600,0)) [1] 0.4531800 0.3559702 0.1613186 > ## Clean up > ## Not run: > ##D rm(gastric,gastric.fit) > ##D dev.off() > ## End(Not run) > > ## THE FOLLOWING EXAMPLE IS TIME-CONSUMING > ## Proportional hazard model > ## Not run: > ##D data(stan) > ##D stan.fit <- sshzd(Surv(futime,status)~futime+age,data=stan) > ##D ## Evaluate fitted hazard > ##D hzdrate.sshzd(stan.fit,data.frame(futime=c(10,20),age=c(20,30))) > ##D ## Plot lambda(t,age=20) > ##D tt <- seq(0,60,leng=101) > ##D hh <- hzdcurve.sshzd(stan.fit,tt,data.frame(age=20)) > ##D plot(tt,hh,type="l") > ##D ## Clean up > ##D rm(stan,stan.fit,tt,hh) > ##D dev.off() > ## End(Not run) > > > > ### *