Comparison of permutation methods

normal permutations

set.seed(42)
perms0<-scanone(fake.f2, pheno.col="phenotype", addcovar=covar,
                intcovar=covar, perm.strata=covar[,1],
                n.perm=100, verbose = F)
summary(perms0)

## LOD thresholds (100 permutations)
##      lod
## 5%  4.10
## 10% 3.89

GWERk scanone permutations

set.seed(42)
perms1<-scanone.GWERk(fake.f2, pheno.col="phenotype", 
                      addcovar=covar, intcovar=covar, 
                      perm.strata=covar[,1],
                      n.perm=100, GWERk=1, verbose = F)
summary(perms1)

## LOD thresholds (100 permutations)
##      lod
## 5%  3.19
## 10% 3.01

plot(as.numeric(perms0), as.numeric(perms1), 
     xlab="standard perms", ylab = "GWER perms", 
     main = "scanone permutation comparison")
abline(a=0,b=1, lty=3)

GWER k=1 permutations

set.seed(42)
perms2.1<-scantwo.GWERk1(fake.f2, chr = c(1,13), 
                         pheno.col="phenotype", 
                         addcovar=covar, intcovar=covar, 
                         perm.strata = covar[,1],
                         n.perm=50, verbose = F)
summary(perms2.1)

## phenotype (50 permutations)
##     full  fv1  int  add  av1  one
## 5%  7.82 5.25 3.38 4.73 2.12 3.32
## 10% 7.11 4.78 3.21 4.46 2.04 2.75

plot(as.numeric(perms2.0$full), as.numeric(perms2.1$full), 
     xlab="standard perms", ylab = "GWER perms", 
     main = "scantwo permutation comparison")
abline(a=0,b=1, lty=3)

Normal stepwiseqtl

pens = calc.penalties(perms2.1, alpha = 0.2)

outsw <- stepwiseqtl(fake.f2, pheno.col="phenotype", chr = c(1,13),
                     max.qtl=3, method="hk", keeptrace=TRUE, covar=covar,
                     penalties = pens, verbose = F)
summary(outsw)

##   QTL object containing genotype probabilities. 
## 
##       name chr   pos n.gen
## Q1  1@37.1   1 37.11     3
## Q2 13@24.0  13 24.03     3
## 
##   Formula: y ~ covar + Q1 + Q2 
## 
##   pLOD:  7.511

summary(fitqtl(fake.f2, pheno.col="phenotype", covar=covar,
               formula=formula(outsw), qtl = outsw, 
               method="hk"))

## 
##      fitqtl summary
## 
## Method: Haley-Knott regression 
## Model:  normal phenotype
## Number of observations : 200 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ covar + Q1 + Q2 
## 
##        df       SS        MS      LOD     %var Pvalue(Chi2)   Pvalue(F)
## Model   5 1168.316 233.66324 12.83204 25.58172 1.870459e-11 3.53888e-11
## Error 194 3398.680  17.51897                                           
## Total 199 4566.996                                                     
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##         df Type III SS    LOD    %var F value Pvalue(Chi2) Pvalue(F)    
## covar    1       34.96 0.4444  0.7654   1.995        0.153     0.159    
## 1@37.1   2      379.29 4.5948  8.3049  10.825        0.000  3.49e-05 ***
## 13@24.0  2      682.57 7.9483 14.9458  19.481        0.000  1.95e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

stepwiseqtl forcing covariate interaction

outsw.int <- stepwiseqtl.int(fake.f2, pheno.col="phenotype", chr = c(1,13),
                             max.qtl=3, method="hk", keeptrace=TRUE, 
                             covar=covar, verbose = F)
summary(outsw.int)

##   QTL object containing genotype probabilities. 
## 
##       name chr   pos n.gen
## Q1  1@37.1   1 37.11     3
## Q2 13@24.0  13 24.03     3
## 
##   Formula: y ~ covar + Q1 + Q2 + Q1:covar + Q2:covar 
## 
##   pLOD:  6.606

summary(fitqtl(fake.f2, pheno.col="phenotype", 
               formula=formula(outsw.int), qtl = outsw.int, 
               covar = covar, method="hk"))

## 
##      fitqtl summary
## 
## Method: Haley-Knott regression 
## Model:  normal phenotype
## Number of observations : 200 
## 
## Full model result
## ----------------------------------  
## Model formula: y ~ covar + Q1 + Q2 + covar:Q1 + covar:Q2 
## 
##        df       SS        MS      LOD     %var Pvalue(Chi2)    Pvalue(F)
## Model   9 1259.581 139.95341 14.01419 27.58007 1.774455e-10 4.601156e-10
## Error 190 3307.416  17.40745                                            
## Total 199 4566.996                                                      
## 
## 
## Drop one QTL at a time ANOVA table: 
## ----------------------------------  
##               df Type III SS     LOD     %var F value Pvalue(Chi2)
## covar          5      126.22 1.62654  2.76373 1.45018        0.187
## 1@37.1         4      465.54 5.71934 10.19365 6.68599        0.000
## 13@24.0        4      684.56 8.16992 14.98930 9.83143        0.000
## covar:1@37.1   2       89.57 1.16049  1.96123 2.57273        0.069
## covar:13@24.0  2        2.92 0.03832  0.06394 0.08387        0.916
##               Pvalue(F)    
## covar             0.208    
## 1@37.1         4.69e-05 ***
## 13@24.0        3.00e-07 ***
## covar:1@37.1      0.079 .  
## covar:13@24.0     0.920    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par(mfrow=c(2,1))
plotLodProfile(outsw, main = "normal stepwise with additive covariate")
plotLodProfile(outsw.int, main = "normal stepwise with interactive covariate")