R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.0 Patched (2005-05-12), 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.

> ### * <HEADER>
> ###
> 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("kinship-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('kinship')
Loading required package: survival
Loading required package: splines
Loading required package: nlme

Attaching package: 'nlme'


	The following object(s) are masked from package:stats :

	 contr.SAS 

Loading required package: lattice
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "as.matrix.bdsmatrix"
> 
> ### * as.matrix.bdsmatrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: as.matrix.bdsmatrix
> ### Title: a function for bdsmatrix
> ### Aliases: as.matrix.bdsmatrix
> ### Keywords: array
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # The function is currently defined as
> ##D function(x)
> ##D {
> ##D         if(class(x) != "bdsmatrix")
> ##D                 stop("argument must be a bdsmatrix object")
> ##D         dd <- dim(x)
> ##D         d3 <- sum(x@blocksize)
> ##D         # dim of square portion
> ##D         d4 <- sum(x@blocksize^2)
> ##D         # size of x@blocks
> ##D         newmat <- matrix(0., dd[1], dd[2], dimnames = x@.Dimnames)
> ##D         temp <- .C("bdsmatrix_index1",
> ##D                 as.integer(length(x@blocksize)),
> ##D                 as.integer(x@blocksize),
> ##D                 as.integer(c(1, 0, 0)),
> ##D                 as.integer(d3),
> ##D                 as.integer(1:d3 - 1),
> ##D                 indexa = integer(d3 * d3),
> ##D                 indexb = 0,
> ##D                 indexc = 0)$indexa
> ##D         newmat[x@permute, x@permute] <- c(x@offdiag, x@blocks)[1 + temp]
> ##D         if(length(x@rmat) > 0) {
> ##D                 newmat[,  - (1:d3)] <- x@rmat
> ##D                 newmat[ - (1:d3),  ] <- t(x@rmat)
> ##D         }
> ##D         newmat
> ##D }
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "bdsBlock"
> 
> ### * bdsBlock
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bdsBlock
> ### Title: Block diagonal matrices.
> ### Aliases: bdsBlock
> ### Keywords: array
> 
> ### ** Examples
> 
> ## Not run: 
> ##D id    <- letters[1:10]
> ##D group <- c(1,1,3,2,3,3,2,3,2,4)
> ##D bdsBlock(id, group)
> ##D     a b d g i c e f h j 
> ##D   a 1 1 0 0 0 0 0 0 0 0
> ##D   b 1 1 0 0 0 0 0 0 0 0
> ##D   d 0 0 1 1 1 0 0 0 0 0
> ##D   g 0 0 1 1 1 0 0 0 0 0
> ##D   i 0 0 1 1 1 0 0 0 0 0
> ##D   c 0 0 0 0 0 1 1 1 1 0
> ##D   e 0 0 0 0 0 1 1 1 1 0
> ##D   f 0 0 0 0 0 1 1 1 1 0
> ##D   h 0 0 0 0 0 1 1 1 1 0
> ##D   j 0 0 0 0 0 0 0 0 0 1
> ##D 
> ##D # Create the matrices for a sparse nested fit of family within city
> ##D group <- paste(mydata$city, mydata$family, sep='/')
> ##D mat1 <- bdsI(group)
> ##D mat2 <- bdsBlock(group, mydata$city)
> ##D fit <- coxme(Surv(time, status) ~ age + sex, data=mydata,
> ##D                random= ~1|group, varlist=list(mat1, mat2))
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "bdsI"
> 
> ### * bdsI
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bdsI
> ### Title: Sparse identity matrices
> ### Aliases: bdsI
> ### Keywords: array
> 
> ### ** Examples
> 
> imat <- bdsI(1:10)
> 
> 
> 
> cleanEx(); ..nameEx <- "bdsmatrix"
> 
> ### * bdsmatrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bdsmatrix
> ### Title: Create a sparse symmetric block diagonal matrix object
> ### Aliases: bdsmatrix print.bdsmatrix
> ### Keywords: manip
> 
> ### ** Examples
> 
> # The matrix shown above is created by
> tmat <- bdsmatrix(c(2,3), c(1,2,1, 3,1,2, 4,3, 5),
+                   rmat=matrix(c(4,6,8,1,2,7,6, 5,7,8,1,2,6,9), ncol=2))
> 
> # Note that only the lower part of the blocks is needed, however, the
> #  entire block set is also allowed, i.e., c(1,2,2,1, 3,1,2,1,4,3,2,3,5)
> 
> 
> 
> cleanEx(); ..nameEx <- "bdsmatrix.ibd"
> 
> ### * bdsmatrix.ibd
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bdsmatrix.ibd
> ### Title: Create a bdsmatrix from a list
> ### Aliases: bdsmatrix.ibd
> ### Keywords: array
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ibdmat <- bdsmatrix.ibd(i,j, ibdval, idlist=subject)
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "besthint"
> 
> ### * besthint
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: besthint
> ### Title: Create a hints matrix for a pedigree.
> ### Aliases: besthint
> ### Keywords: array
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # Find a good plot, only trying to avoid dotted connectors
> ##D myped$hints <- besthint(myped, wt=c(1000,100,0))
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "coxme"
> 
> ### * coxme
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: coxme
> ### Title: Fit a mixed-effects Cox model
> ### Aliases: coxme coxme.fit coxme.varcheck print.coxme
> ### Keywords: survival
> 
> ### ** Examples
> 
> ## Not run: 
> ##D coxme(Surv(time, status) ~ rx, data=rats, random= ~1|litter)
> ##D 
> ##D Cox mixed-effects kinship model fit by maximum likelihood
> ##D   Data: rats 
> ##D   n= 150 
> ##D                     NULL Integrated Penalized 
> ##D Log-likelihood -185.6556   -180.849  -173.774
> ##D 
> ##D   Penalized loglik: chisq= 23.76 on 13.17 degrees of freedom, p= 0.036 
> ##D  Integrated loglik: chisq= 9.61 on 2 degrees of freedom, p= 0.0082 
> ##D 
> ##D Fixed effects: Surv(time, status) ~ rx 
> ##D         coef exp(coef)  se(coef)        z         p 
> ##D rx 0.9132825  2.492491 0.3226856 2.830255 0.0046511
> ##D 
> ##D Random effects:  ~ 1 | litter 
> ##D              litter 
> ##D Variance: 0.4255484
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "familycheck"
> 
> ### * familycheck
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: familycheck
> ### Title: Error check for a family classification
> ### Aliases: familycheck
> ### Keywords: manip
> 
> ### ** Examples
> 
> ## Not run: 
> ##D #
> ##D # This is from a pedigree that had some identifier errors
> ##D #
> ##D > checkit<- familycheck(ids2$famid, ids2$gid, ids2$fatherid, ids2$motherid)
> ##D > table(checkit$split)    # should be all 1's
> ##D    0   1 2 
> ##D  112 424 4
> ##D # Shows 112 of the "families" were actually isolated individuals,
> ##D #  and that 4 of the families actually split into 2.  
> ##D #  In one case, a mistyped father id caused one child, along with his spouse
> ##D #  and children, to be "set adrift" from the connected pedigree.
> ##D 
> ##D > table(checkit$join)
> ##D    0 1 2 
> ##D  531 6 3
> ##D #
> ##D # There are 6 families with 1 other joined to them (3 pairs), and 3 with
> ##D #  2 others added to them (one triplet).
> ##D # For instance, a single mistyped father id of someone in family 319,
> ##D #  which was by bad luck the id of someone else in family 339,
> ##D #  was sufficient to join two groups.
> ##D > attr(checkit, 'join')
> ##D     [,1] [,2] [,3] [,4] [,5] [,6] [,7] 
> ##D  31   78    0    0    0    0    0    0
> ##D  32    3   15    0    0    0    0    0
> ##D  33    6    0   12    0    0    0    0
> ##D  63    0    0    0   63    0    0    0
> ##D  65    0    0    0   17   16    0    0
> ##D 122    0    0    0    0    0   16    0
> ##D 127    0    0    0    0    0   30    0
> ##D 319    0    0    0    0    0    0   20
> ##D 339    0    0    0    0    0    0   37
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "gchol"
> 
> ### * gchol
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: gchol
> ### Title: Generalized Cholesky decompostion
> ### Aliases: gchol gchol.bdsmatrix
> ### Keywords: array
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # Create a matrix that is symmetric, but not positive definite
> ##D #   The matrix temp has column 6 redundant with cols 1-5
> ##D smat <- matrix(1:64, ncol=8)
> ##D smat <- smat + t(smat) + diag(rep(20,8))  #smat is 8 by 8 symmetric
> ##D temp <-  smat[c(1:5, 5:8), c(1:5, 5:8)]
> ##D ch1  <- gchol(temp)
> ##D 
> ##D print(as.matrix(ch1))   # print out L
> ##D print(diag(ch1))        # print out D
> ##D aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
> ##D aeq(diag(ch1)[6], 0)    # Check that it has a zero in the proper place
> ##D 
> ##D ginv <- solve(ch1)    # see if I get a generalized inverse
> ##D aeq(temp %*% ginv %*% temp, temp)
> ##D aeq(ginv %*% temp %*% ginv, ginv)
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "kinship"
> 
> ### * kinship
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: kinship
> ### Title: Compute a kinship matrix
> ### Aliases: kinship kindepth
> ### Keywords: array
> 
> ### ** Examples
> 
> ## Not run: 
> ##D test1 <- data.frame(id  =c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14),
> ##D                     mom =c(0, 0, 0, 0, 2, 2, 4, 4, 6,  2,  0,  0, 12, 13),
> ##D                     dad =c(0, 0, 0, 0, 1, 1, 3, 3, 3,  7,  0,  0, 11, 10),
> ##D                     sex =c(0, 1, 0, 1, 0, 1, 0, 1, 0,  0,  0,  1,  1,  1))
> ##D round(8*kinship(test1$id, test1$dad, test1$mom))
> ##D 
> ##D    1 2 3 4 5 6 7 8 9 10 11 12 13 14 
> ##D  1 4 0 0 0 2 2 0 0 1  0  0  0  0  0
> ##D  2 0 4 0 0 2 2 0 0 1  2  0  0  0  1
> ##D  3 0 0 4 0 0 0 2 2 2  1  0  0  0  0
> ##D  4 0 0 0 4 0 0 2 2 0  1  0  0  0  0
> ##D  5 2 2 0 0 4 2 0 0 1  1  0  0  0  0
> ##D  6 2 2 0 0 2 4 0 0 2  1  0  0  0  0
> ##D  7 0 0 2 2 0 0 4 2 1  2  0  0  0  1
> ##D  8 0 0 2 2 0 0 2 4 1  1  0  0  0  0
> ##D  9 1 1 2 0 1 2 1 1 4  1  0  0  0  0
> ##D 10 0 2 1 1 1 1 2 1 1  4  0  0  0  2
> ##D 11 0 0 0 0 0 0 0 0 0  0  4  0  2  1
> ##D 12 0 0 0 0 0 0 0 0 0  0  0  4  2  1
> ##D 13 0 0 0 0 0 0 0 0 0  0  2  2  4  2
> ##D 14 0 1 0 0 0 0 1 0 0  2  1  1  2  4
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "lmekin"
> 
> ### * lmekin
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: lmekin
> ### Title: Linear Mixed Effects model using a kinship matrix.
> ### Aliases: lmekin print.lmekin getCovariateFormula2 getCrossedTerms
> ###   getGroupsFormula2
> ### Keywords: regression
> 
> ### ** Examples
> 
> ## Not run: 
> ##D #
> ##D # Make a kinship matrix for the entire study
> ##D # These two functions are NOT fast, the makekinship one in particular
> ##D #
> ##D cfam <- makefamid(main$gid, main$momid, main$dadid)
> ##D kmat <- makekinship(cfam, main$gid, main$momid, main$dadid)
> ##D 
> ##D # The kinship matrix for the females only: quite a bit smaller
> ##D #
> ##D kid <- dimnames(kmat)[[1]]
> ##D temp <- main$sex[match(kid, main$gid)] == 'F'
> ##D fkmat <- kmat[temp,temp]
> ##D 
> ##D # The dimnames on kmat are the gid value, which are necessary to match
> ##D #  the appropriate row/col of kmat to the analysis data set
> ##D # A look at %dense tissue on a mammogram, with age at mammogram and
> ##D #  weight as covariates, and a familial random effect
> ##D #
> ##D fit <- lmekin(percdens ~ mammage + weight, data=anal1,
> ##D                  random = ~1|gid, varlist=list(fkmat))
> ##D 
> ##D Linear mixed-effects kinship model fit by maximum likelihood
> ##D   Data: anal1 
> ##D   Log-likelihood = -6093.917 
> ##D   n= 1535 
> ##D 
> ##D Fixed effects: percdens ~ mammage + weight 
> ##D  (Intercept)    mammage     weight 
> ##D      87.1593 -0.5333198 -0.1948871
> ##D 
> ##D Random effects:  ~ 1 | gid 
> ##D          Kinship Residual 
> ##D StdDev: 7.801603 10.26612
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "makefamid"
> 
> ### * makefamid
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: makefamid
> ### Title: Identify family groups
> ### Aliases: makefamid
> ### Keywords: manip
> 
> ### ** Examples
> 
> ## Not run: 
> ##D > newid <- makefamid(cdata$gid, cdata$dadid, cdata$momid)
> ##D > table(newid==0)
> ##D  FALSE TRUE 
> ##D  17859 8191
> ##D # So nearly 1/3 of the individuals are not blood relatives.
> ##D 
> ##D > kin1 <- makekinship(cdata$famid, cdata$gid, cdata$dadid, cdata$momid)
> ##D > kin2 <- makekinship(newid, cdata$gid, cdata$dadid, cdata$momid, unique=0)
> ##D > dim(kin2)
> ##D [1] 26050 26050
> ##D > dim(kin1)
> ##D [1] 26050 26050
> ##D 
> ##D > length(kin2@blocks)/length(kin1@blocks)
> ##D [1] 0.542462
> ##D # Basing kin1 on newid rather than cdata$famid (where marry-ins were each
> ##D #   labeled as members of one of the 426 families) reduced its size by just
> ##D #   less than half.
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "makekinship"
> 
> ### * makekinship
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: makekinship
> ### Title: Create a sparse kinship matrix
> ### Aliases: makekinship
> ### Keywords: manip
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # Data set from a large family study of breast cancer
> ##D #  there are 26050 subjects in the file, from 426 families
> ##D > table(cdata$sex)
> ##D      F     M 
> ##D  12699 13351
> ##D > length(unique(cdata$famid))
> ##D [1] 426
> ##D 
> ##D > kin1 <- makekinship(cdata$famid, cdata$gid, cdata$dadid, cdata$momid)
> ##D > dim(kin1)
> ##D [1] 26050 26050
> ##D > class(kin1)
> ##D [1] "bdsmatrix"
> ##D # The next line shows that few of the elements of the full matrix are >0
> ##D > length(kin1@blocks)/ prod(dim(kin1))
> ##D [1] 0.00164925
> ##D 
> ##D # kinship matrix for the females only
> ##D > femid <- cdata$gid[cdata$sex=='F']
> ##D > femindex <- !is.na(match(dimnames(kin1)[[1]], femid))
> ##D > kin2 <- kin1[femindex, femindex]
> ##D #
> ##D # Note that "femindex <- match(femid, dimnames(kin1)[[1]])" is wrong, since
> ##D #  then kin1[femindex, femindex] might improperly reorder the rows/cols 
> ##D #  (if families were not contiguous in cdata).  
> ##D # However sort(match(femid, dimnames(kin1)[[1]])) would be okay.
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "pedigree"
> 
> ### * pedigree
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: pedigree
> ### Title: Create pedigree structure
> ### Aliases: pedigree
> ### Keywords: dplot
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ptemp <- pedigree(id=d10$upn, dadid=d10$dadid,momid=d10$momid, 
> ##D                   sex=d10$sex, affected=d10$affect)
> ##D plot(ptemp)
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.pedigree"
> 
> ### * plot.pedigree
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.pedigree
> ### Title: plot pedigrees
> ### Aliases: plot.pedigree
> ### Keywords: hplot
> 
> ### ** Examples
> 
> ## Not run: 
> ##D ptemp <- pedigree(id=d10$upn, dadid=d10$dadid, momid=d10$momid, 
> ##D                   sex=d10$sex, affected=d10$affect)
> ##D 
> ##D plot(ptemp)
> ##D 
> ##D col.founder <- rep(1,length(d10$affect))
> ##D col.founder[d10$id==1] <- 2
> ##D 
> ##D plot(ptemp, affected=cbind(d10$affect,d10$affect2,d10$affect3,d10$affect4),
> ##D      col=col.founder, id=paste(ptemp$id,'\n','(',d10$post,')',sep=''),
> ##D      angle=c(90,80,70,60), density=c(-1,90,70,50))
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "solve.bdsmatrix"
> 
> ### * solve.bdsmatrix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: solve.bdsmatrix
> ### Title: Solve a matrix equation using the generalized Cholesky
> ###   decompostion
> ### Aliases: solve.bdsmatrix
> ### Keywords: algebra
> 
> ### ** Examples
> 
> ## Not run: 
> ##D tmat <- bdsmatrix(c(3,2,2,4), 
> ##D               c(22,1,2,21,3,20,19,4,18,17,5,16,15,6,7, 8,14,9,10,13,11,12),
> ##D               matrix(c(1,0,1,1,0,0,1,1,0,1,0,10,0,
> ##D                        0,1,1,0,1,1,0,1,1,0,1,0,10), ncol=2))
> ##D dim(tmat)
> ##D solve(tmat, cbind(1:13, rep(1,13)))
> ## End(Not run)
> 
> 
> 
> cleanEx(); ..nameEx <- "solve.gchol"
> 
> ### * solve.gchol
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: solve.gchol
> ### Title: Solve a matrix equation using the generalized Cholesky
> ###   decompostion
> ### Aliases: solve.gchol
> ### Keywords: algebra
> 
> ### ** Examples
> 
> ## Not run: 
> ##D # Create a matrix that is symmetric, but not positive definite
> ##D #   The matrix temp has column 6 redundant with cols 1-5
> ##D smat <- matrix(1:64, ncol=8)
> ##D smat <- smat + t(smat) + diag(rep(20,8))  #smat is 8 by 8 symmetric
> ##D temp <-  smat[c(1:5, 5:8), c(1:5, 5:8)]
> ##D ch1  <- gchol(temp)
> ##D 
> ##D print(as.matrix(ch1))   # print out L
> ##D print(diag(ch1))        # print out D
> ##D aeq <- function(x,y) all.equal(as.vector(x), as.vector(y))
> ##D aeq(diag(ch1)[6], 0)    # Check that it has a zero in the proper place
> ##D 
> ##D ginv <- solve(ch1)    # see if I get a generalized inverse
> ##D aeq(temp %*% ginv %*% temp, temp)
> ##D aeq(ginv %*% temp %*% ginv, ginv)
> ## End(Not run)
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  2.49 0.09 2.64 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')