Nothing
f.design.make <- function(maternal, response, info, ret.characteristics = F){
#
# THE PROGRAM ESTIMATES EFFECTS OF SEVERAL ALLELES IN A CASE-TRIAD, CASE-CONTROL-TRIAD OR CASE-CONTROL DESIGN
# CREATES AN APPROPRIATE DESIGN MATRIX FOR USE IN f.tri.glm
#
# THE MODEL ASSUMES HARDY-WEINBERG EQUILIBRIUM (AND RARE DISEASE, IF NECESSARY, IN CC)
# IMPORTANT: OBSERVED FREQUENCIES (observed) MUST COME FROM A COMPLETE
# GRID IN APPROPRIATE ORDER! THE ORDER IS:
# TRIAD: BY MOTHER'S FIRST ALLELE AND SECOND ALLELE,
# THEN THE FATHER'S FIRST AND SECOND
# CASE-CONTROL-TRIAD: BY MOTHER'S FIRST ALLELE AND SECOND ALLELE,
# THEN THE FATHER'S FIRST AND SECOND, AND FINALLY CONTROL (1) AND CASE (2),
# SO THAT TOP HALF OF THE DATA ARE CONTROLS
# CASE-CONTROL: BY FIRST ALLELE, SECOND ALLELE
# AND FINALLY CONTROL (1) AND CASE (2), SO THAT TOP HALF OF THE DATA ARE
# CONTROLS
#
# maternal DETERMINES WHETHER OR NOT MATERNAL EFFECTS SHOULD BE ESTIMATED
# (NOT AVAILABLE FOR THE CASE-CONTROL DESIGN)
#
#
## EXTRACT VARIABLES
.ref.cat <- info$haplos$ref.cat
.xchrom <- info$model$xchrom
.design <- info$model$design
.poo <- info$model$poo
.hwe <- info$model$hwe
if(is.null(.hwe)) .hwe <- TRUE
.sel.sex <- info$variables$sel.sex
.sel.sex.yes <- !is.null(.sel.sex)
.comb.sex <- info$model$comb.sex
.covar.yes <- !is.null(info$variables$covar)
.A <- sum(info$haplos$selected.haplotypes) ## NUMBER OF ALLELES (HAPLOTYPES)
#
###########
# CREATE "BASE" GRID OF PARENTAL GENOTYPES (AND POSSIBLY SEX AND/OR CASE-CONTROL STATUS)
###########
.tmp <- factor(1:.A)
.mf.tmp <- list(m1 = .tmp, m2 = .tmp, f1 = .tmp, f2 = .tmp, sex = c(0,1), cc = c(0,1))
## HERE, 0 AND 1 (MALES AND FEMALES, RESP.) ARE USED AS DUMMY FOR sex. ORIGINAL CODING (OUTSIDE f.tri.glm AND f.make.design) SHOULD BE 1 (MALES) AND 2 (FEMALES)
## HERE, 0 AND 1 (CONTROL AND CASE, RESP.) ARE USED AS DUMMY FOR cc. ORIGINAL CODING (OUTSIDE f.tri.glm AND f.make.design) SHOULD BE 1 (CONTROL) AND 2 (CASE). IN INPUT FILE, THE LARGEST IS ALWAYS THE CASE
.cond.triad <- (.design %in% c("triad", "cc.triad"))
.cond.f1 <- !.xchrom & (.design != "cc")
.cond.sex <- .xchrom & !.sel.sex.yes
.cond.cc <- (.design %in% c("cc", "cc.triad"))
#
## m1
if(!.cond.triad) .mf.tmp[["m1"]] <- NULL
## m2
## f1
if(!.cond.f1) .mf.tmp[["f1"]] <- NULL
## f2
## sex
if(!.cond.sex) .mf.tmp[["sex"]] <- NULL
## cc
if(!.cond.cc) .mf.tmp[["cc"]] <- NULL
#
##
#if((.design == "cc") & .xchrom) stop("Not implemented")## MERK: DE OVENSTAAENDE SELEKSJONENE SKAL FUNGERE FOR DENNE SITUASJONEN OGSAA!
if(.poo & (.design == "cc"))stop("Not implemented")## DETTE BLIR DET OGSAA SJEKKET FOR I f.check.pars
if(.covar.yes) .mf.tmp$covar <- factor(seq(along = info$variables$covar.codes))
#
## EXPAND GRID
.mf <- do.call("expand.grid", .mf.tmp)
if(.xchrom){
# Define a sex variable for xchrom
#
## RECALL THAT sex IS NOT PART OF .mf IF ONLY ONE SEX IS SELECTED
if(.comb.sex %in% c("single", "double")){
.sex <- .mf$sex
}
if(.comb.sex == "males"){
.sex <- rep(0, nrow(.mf))
}
if(.comb.sex == "females"){
.sex <- rep(1, nrow(.mf))
}
}
if(.xchrom & (.design == "cc")){
## For this special case, reduce grid
## MERK: dette kunne ogsaa vaert gjort for alle xchrom, men de er i data lagt ut som ambiguities
.rem <- (.sex == 0) & (.mf$m2 != .mf$f2)
.mf <- .mf[!.rem,]
.sex <- .sex[!.rem]
}
if(ret.characteristics){
## Return what's needed to regenerate grid
if(.xchrom & (.design == "cc")){
## For the time being, return entire exapnded and reduced grid for matching. Remember to change m2, f2 into c1, c2, and sex <- sex + 1
return(.mf)
}else{
.char <- sapply(.mf.tmp, length)
return(.char)
}
}
#
###########
# CREATE DUMMY VARIABLES FOR PARENTAL GENOTYPES
###########
if(.cond.triad){
.m1.dum <- model.matrix( ~ -1 + m1, data = .mf)
}
.m2.dum <- model.matrix( ~ -1 + m2, data = .mf)
if(.cond.f1){
.f1.dum <- model.matrix( ~ -1 + f1, data = .mf)
}
.f2.dum <- model.matrix( ~ -1 + f2, data = .mf)
#
###########
# CREATE COUNTING VARIABLES FOR PARENTS AND CHILDREN
###########
## MOTHER
if(.cond.triad){
.m.dumsum <- .m1.dum + .m2.dum
dimnames(.m.dumsum)[[2]] <- paste("m", 1:.A, sep = "")
}else{
.m.dumsum <- .m2.dum
}
## FATHER
if(.cond.f1){
.f.dumsum <- .f1.dum + .f2.dum
dimnames(.f.dumsum)[[2]] <- paste("f", 1:.A, sep = "")
}else{
.f.dumsum <- .f2.dum
}
#
## PARENTS COMBINED, FOR ESTIMATING ALLELE FREQUENCIES
.parents.dumsum <- .m.dumsum + .f.dumsum
dimnames(.parents.dumsum)[[2]] <- paste("mf", 1:.A, sep = "")
#
## CHILD, DEFINE .c.dumsum:
if(.xchrom){
## X-INACTIVATION(comb.sex == "double"): BOYS GET TWICE THE EFFECT. GIRLS GET ONE OR TWO.
## MEANS THAT SINGLE DOSE IN BOYS CORRESPONDS TO DOUBLE DOSE IN GIRLS,
## SINGLE DOSE IN GIRLS IS SQUARE ROOT OF DOUBLE DOSE (WHEN response = "mult")
## BOYS: 1, RR^2, GIRLS: 1, RR, RR^2 (response = "mult")
## BOYS: 1, RR1^2*RR2, GIRLS: 1, RR1, RR1^2*RR2 (response = "free")
#
## DOSE-RESPONSE(comb.sex == "single"): BOYS GET SINGLE DOSE EFFECT. GIRLS GET ONE OR TWO.
## MEANS THAT SINGLE DOSE IN BOYS CORRESPONDS TO SINGLE DOSE IN GIRLS,
## DOUBLE DOSE IN GIRLS IS SQUARE OF SINGLE DOSE. (WHEN response = "mult")
## BOYS: 1, RR, GIRLS: 1, RR, RR^2 (response = "mult")
## BOYS: 1, RR1, GIRLS: 1, RR1, RR1^2*RR2 (response = "free")
#
.m2.girls.dum <- .m2.dum * .sex # ALLELE FROM MOTHER, IN GIRLS
.f2.girls.dum <- .f2.dum * .sex # ALLELE FROM FATHER, IN GIRLS
.m2.boys.dum <- .m2.dum * (1 - .sex) # ALLELE FROM MOTHER, IN BOYS
#
## DEFINE COMBINATION
if(.comb.sex == "double"){
.c.dumsum <- .m2.girls.dum + .f2.girls.dum + 2 * .m2.boys.dum
}else{
.c.dumsum <- .m2.girls.dum + .f2.girls.dum + 1 * .m2.boys.dum
}
if(.design == "cc"){
## Use only maternal allele for frequency computation
.parents.dumsum[.sex == 0, ] <- .m.dumsum[.sex == 0, ]
}
}else{
## CHILD, STANDARD AUTOSOMAL MODEL
.c.dumsum <- .m2.dum + .f2.dum
}
dimnames(.c.dumsum)[[2]] <- paste("c", 1:.A, sep = "")
#
## WHEN RESPONSE IS FREE, ADD STAR TERM
if(response == "free"){
## SPECIAL CODING FOR HOMOZYGOTES:
.c.dd <- (.c.dumsum == 2) + 0 # MODEL: R^2*Rstar
# SET APPROPRIATE COLUMN NAMES:
dimnames(.c.dd)[[2]] <- paste("cdd", 1:.A, sep = "") #
#
if(maternal){
.m.dd <- (.m.dumsum == 2) + 0 # MODEL: R^2*Rstar FOR MAT. EFF
# SET APPROPRIATE COLUMN NAMES:
dimnames(.m.dd)[[2]] <- paste("mdd", 1:.A, sep = "")
}# END if(maternal)
}# END if(response == "free")
#
##
if(.covar.yes){
.tmpmat <- cbind(.parents.dumsum, .mf[, c("covar"), drop = F])
.form <- paste("mf", 1:.A, sep = "", collapse = " + ")
.form <- paste("~ -1 + (", .form, "):covar", sep = "")
.form <- formula(.form)
.parents.dumsum <- model.matrix(.form, data = .tmpmat)
cat("kontroller at sortering etc. blir rett!\n")
.navn <- colnames(.parents.dumsum)
.navn <- gsub(":", "_", .navn)
# .navn <- gsub("_covar1", "", .navn)
# cat("ad hoc!\n")
colnames(.parents.dumsum) <- .navn
}
# SET UP FINAL DESIGN MATRIX, REMOVING REDUNDANT .ref.cat-ROW.
# NOTE: DIALLELIC SITUATION REQUIRES THE REMOVAL OF ONE DOUBLE DOSE PARAMETER:
#
## START WITH sex AND cc (WHICH ARE NULL IF NOT RELEVANT)
.design.matrix0 <- cbind(sex = .mf$sex, cc = .mf$cc)
## ADD .parents.dumsum, FOR ESTIMATING HAPLOTYPE FREQUENCIES
.design.matrix0 <- cbind(.design.matrix0, .parents.dumsum)
## KEEP THE EFFECTS PART SEPARATE TO BEGIN WITH, IN CASE DESIGN INCLUDES cc
.design.matrix1 <- NULL
#
## IF response IS MORE THAN simple, ADD THE NECESSARY DESIGN VARIABLES
if(response != "simple"){
if(.poo){
if(.xchrom){
if(.comb.sex == "single"){
## A (SINGLE) DOSE IN MALES CORRESPONDS TO A SINGLE, MATERNALLY INHERITED, DOSE IN FEMALES
.mat <- .m2.girls.dum + 1 * .m2.boys.dum
.pat <- .f2.girls.dum
}
if(.comb.sex == "double"){
## A (SINGLE) DOSE IN MALES CORRESPONDS TO A DOUBLE DOSE IN FEMALES, I.E. A SINGLE MATERNALLY INHERITED + A SINGLE PATERNALLY INHERITED + THE Rstar DOUBLE DOSE
## SO THAT A (SINGLE) DOSE IN MALES CORRESPONDS TO THE "FULL PACKAGE" IN FEMALES
.mat <- .m2.girls.dum + 1 * .m2.boys.dum
.pat <- .f2.girls.dum + 1 * .m2.boys.dum
}
if(.comb.sex == "males"){
.mat <- .m2.boys.dum
.pat <- NULL
}
if(.comb.sex == "females"){
.mat <- .m2.girls.dum
.pat <- .f2.girls.dum
}
## RENAME COLUMNS AND DROP REF CAT
colnames(.mat) <- sub("m2", "cm", colnames(.mat))
.mat <- .mat[, - .ref.cat, drop = F]
if(.comb.sex != "males"){
colnames(.pat) <- sub("f2", "cf", colnames(.pat))
.pat <- .pat[, - .ref.cat, drop = F]
}
.design.matrix1 <- cbind(.design.matrix1, .mat, .pat)
}else{
## RENAME COLUMNS
colnames(.m2.dum) <- sub("m2", "cm", colnames(.m2.dum))
colnames(.f2.dum) <- sub("f2", "cf", colnames(.f2.dum))
## DROP REF CAT
.m2.dum <- .m2.dum[, - .ref.cat, drop = F]
.f2.dum <- .f2.dum[, - .ref.cat, drop = F]
.design.matrix1 <- cbind(.design.matrix1, .m2.dum, .f2.dum)
}
}else{# if(!.poo)
.c.dumsum <- .c.dumsum[, - .ref.cat, drop = F]
.design.matrix1 <- cbind(.design.matrix1, .c.dumsum)
}
if(response == "free"){
if(.A == 2) .c.dd <- .c.dd[, - .ref.cat, drop = F]
.design.matrix1 <- cbind(.design.matrix1, .c.dd)
}
if(maternal){
.m.dumsum <- .m.dumsum[, - .ref.cat, drop = F]
.design.matrix1 <- cbind(.design.matrix1, .m.dumsum)
if(response == "free"){
if(.A == 2) .m.dd <- .m.dd[, - .ref.cat, drop = F]
.design.matrix1 <- cbind(.design.matrix1, .m.dd)
}
}
}
#
## LET THE FORMULA PART APPLY ONLY FOR CASES.
## (WARNING! RARE DISEASE ASSUMPTION HERE!)
if(.cond.cc){
.design.matrix1 <- .design.matrix1 * .mf$cc
}
#
## JOIN AND CONVERT TO DATA FRAME
.design.matrix <- as.dframe(cbind(.design.matrix0, .design.matrix1))
#
##
# return(cbind(.mf, .design.matrix))
return(.design.matrix)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.