#simulation of effects for cis trans test
#wts<-vector of weights indicating the effect size of the factors in each group
#nind<-number of individuals / group
#dis<-numeric indicator of dispersion around mean
#xbar<-intercept
sim.ct<-function(xbar=5, nind=80, sd.sim=.4,
factors=list(allele=c("hal","fil"),generation=c("F0","F1"),treatment=c("wet","dry")),
model= "y~cis+trans.f+trans.h+trt+cis*trt + trans.f*trt + trans.h*trt",
wts=c(1,0,0,0,0,0,0), type="nbinom"){
terms<-gsub(" ","",model); terms<-gsub("y~","",terms)
terms<-unlist(strsplit(terms,"[+]"))
model<-as.formula(model)
grps<-expand.grid(factors)
rownames(grps)<-paste(grps$allele, grps$generation, grps$treatment, sep=".")
grps$effect<-0
wts.fil.F0.wet<-c(1,1,0,1,1,1,0)
wts.fil.F0.dry<-c(1,1,0,0,0,0,0)
wts.hal.F0.wet<-c(0,0,1,1,0,0,1)
wts.hal.F0.dry<-c(0,0,1,0,0,0,0)
wts.fil.F1.wet<-c(1,1,1,1,1,1,1)
wts.fil.F1.dry<-c(1,1,1,0,0,0,0)
wts.hal.F1.wet<-c(0,1,1,1,0,1,1)
wts.hal.F1.dry<-c(0,1,1,0,0,0,0)
grps["fil.F0.wet","effect"]<-sum(wts*wts.fil.F0.wet)
grps["hal.F0.wet","effect"]<-sum(wts*wts.hal.F0.wet)
grps["fil.F1.wet","effect"]<-sum(wts*wts.fil.F1.wet)
grps["hal.F1.wet","effect"]<-sum(wts*wts.hal.F1.wet)
grps["fil.F0.dry","effect"]<-sum(wts*wts.fil.F0.dry)
grps["hal.F0.dry","effect"]<-sum(wts*wts.hal.F0.dry)
grps["fil.F1.dry","effect"]<-sum(wts*wts.fil.F1.dry)
grps["hal.F1.dry","effect"]<-sum(wts*wts.hal.F1.dry)
grps$means<-xbar+grps$effect
ind.pergroup<-nind/dim(grps)[1]
if(type=="normal"){
for(i in 1:ind.pergroup){
dat<-data.frame(rnorm(8,mean=grps$means,sd=sd.sim))
colnames(dat)<-paste("ind",i,sep=".")
grps<-cbind(grps,dat)
}
}else{
for(i in 1:ind.pergroup){
dat<-data.frame(rnegbin(8,mu=exp(grps$means),theta=exp(sd.sim)))
colnames(dat)<-paste("ind",i,sep=".")
grps<-cbind(grps,dat)
}
}
library(reshape2)
library(ggplot2)
library(MASS)
# Specify id.vars: the variables to keep but not split apart on
grps.long<-melt(grps, id.vars=c("allele","generation","treatment","means","effect"))
grps.long$cis<-grps.long$allele
grps.long$trans.f<-ifelse(grps.long$generation == "F1" | grps.long$allele =="fil","trans.fil","not")
grps.long$trans.h<-ifelse(grps.long$generation == "F1" | grps.long$allele =="hal","trans.hal","not")
grps.long$treatment<-as.factor(grps.long$treatment)
grps.long$trans.h<-as.factor(grps.long$trans.h)
grps.long$trans.f<-as.factor(grps.long$trans.f)
if(type=="normal"){
glm.out<-lm(value~1+cis+trans.f+trans.h+cis*treatment+trans.f*treatment+trans.h*treatment, data=grps.long)
print(summary(glm.out))
print(ggplot(grps.long, aes(x=allele,y=value,color=generation))+
geom_boxplot()+ facet_grid(treatment~.)+theme_bw()
)
}else{
glm.out<-glm.nb(value~cis+trans.f+trans.h+cis*treatment+trans.f*treatment+trans.h*treatment, link="log", data=grps.long)
print(summary(glm.out)$coefficients)
print(ggplot(grps.long, aes(x=allele,y=value,color=generation))+
geom_boxplot()+ facet_grid(treatment~.)+theme_bw()
)
}
return(glm.out)
}
#cis effect only
test<-sim.ct(wts=c(2,0,0,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8881 -0.3353 0.0177 0.2727 1.0520
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.88005 0.22837 21.369 <2e-16 ***
## cisfil 2.25957 0.18647 12.118 <2e-16 ***
## trans.ftrans.fil -0.11458 0.18647 -0.614 0.541
## trans.htrans.hal 0.10461 0.18647 0.561 0.577
## treatmentdry 0.20922 0.32297 0.648 0.519
## cisfil:treatmentdry -0.07559 0.26370 -0.287 0.775
## trans.ftrans.fil:treatmentdry -0.04785 0.26370 -0.181 0.857
## trans.htrans.hal:treatmentdry -0.10229 0.26370 -0.388 0.699
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.417 on 72 degrees of freedom
## Multiple R-squared: 0.8787, Adjusted R-squared: 0.867
## F-statistic: 74.54 on 7 and 72 DF, p-value: < 2.2e-16
#cis and trans.f effect
test<-sim.ct(wts=c(2,2,0,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.88470 -0.27273 0.00328 0.23359 0.99191
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.889050 0.217489 22.480 <2e-16 ***
## cisfil 2.055413 0.177579 11.575 <2e-16 ***
## trans.ftrans.fil 2.118619 0.177579 11.931 <2e-16 ***
## trans.htrans.hal 0.018007 0.177579 0.101 0.920
## treatmentdry 0.125537 0.307575 0.408 0.684
## cisfil:treatmentdry 0.009347 0.251134 0.037 0.970
## trans.ftrans.fil:treatmentdry -0.054419 0.251134 -0.217 0.829
## trans.htrans.hal:treatmentdry -0.245767 0.251134 -0.979 0.331
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3971 on 72 degrees of freedom
## Multiple R-squared: 0.9555, Adjusted R-squared: 0.9512
## F-statistic: 220.8 on 7 and 72 DF, p-value: < 2.2e-16
#cis and trans.h effect
test<-sim.ct(wts=c(2,0,2,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.83754 -0.25966 -0.00887 0.23306 0.94756
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.14856 0.20479 25.140 <2e-16 ***
## cisfil 1.85730 0.16721 11.107 <2e-16 ***
## trans.ftrans.fil 0.01017 0.16721 0.061 0.952
## trans.htrans.hal 2.04152 0.16721 12.209 <2e-16 ***
## treatmentdry -0.05696 0.28962 -0.197 0.845
## cisfil:treatmentdry 0.12740 0.23647 0.539 0.592
## trans.ftrans.fil:treatmentdry -0.17371 0.23647 -0.735 0.465
## trans.htrans.hal:treatmentdry 0.08316 0.23647 0.352 0.726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3739 on 72 degrees of freedom
## Multiple R-squared: 0.8516, Adjusted R-squared: 0.8372
## F-statistic: 59.03 on 7 and 72 DF, p-value: < 2.2e-16
#cis and cis*trt int
test<-sim.ct(wts=c(1,0,0,0,2,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.82095 -0.21736 -0.04014 0.18586 1.20587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.13040 0.23288 22.030 < 2e-16 ***
## cisfil 3.07823 0.19015 16.188 < 2e-16 ***
## trans.ftrans.fil 0.02216 0.19015 0.117 0.908
## trans.htrans.hal -0.26194 0.19015 -1.378 0.173
## treatmentdry -0.09008 0.32935 -0.274 0.785
## cisfil:treatmentdry -2.22630 0.26891 -8.279 4.62e-12 ***
## trans.ftrans.fil:treatmentdry 0.15574 0.26891 0.579 0.564
## trans.htrans.hal:treatmentdry 0.17284 0.26891 0.643 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4252 on 72 degrees of freedom
## Multiple R-squared: 0.9104, Adjusted R-squared: 0.9017
## F-statistic: 104.5 on 7 and 72 DF, p-value: < 2.2e-16
#trans.f only
test<-sim.ct(wts=c(0,2,0,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.35265 -0.18423 -0.03143 0.20818 0.88338
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.95159 0.20211 24.500 <2e-16 ***
## cisfil -0.07422 0.16502 -0.450 0.654
## trans.ftrans.fil 2.12158 0.16502 12.856 <2e-16 ***
## trans.htrans.hal -0.02722 0.16502 -0.165 0.869
## treatmentdry -0.16530 0.28583 -0.578 0.565
## cisfil:treatmentdry 0.05209 0.23338 0.223 0.824
## trans.ftrans.fil:treatmentdry 0.23753 0.23338 1.018 0.312
## trans.htrans.hal:treatmentdry -0.08219 0.23338 -0.352 0.726
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.369 on 72 degrees of freedom
## Multiple R-squared: 0.8845, Adjusted R-squared: 0.8733
## F-statistic: 78.8 on 7 and 72 DF, p-value: < 2.2e-16
#trans.h only
test<-sim.ct(wts=c(0,0,2,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.12654 -0.22660 0.08034 0.28197 0.73119
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.05428 0.21054 24.007 <2e-16 ***
## cisfil -0.15812 0.17190 -0.920 0.3607
## trans.ftrans.fil 0.18533 0.17190 1.078 0.2846
## trans.htrans.hal 1.85518 0.17190 10.792 <2e-16 ***
## treatmentdry -0.06922 0.29774 -0.232 0.8168
## cisfil:treatmentdry 0.45266 0.24311 1.862 0.0667 .
## trans.ftrans.fil:treatmentdry -0.49306 0.24311 -2.028 0.0462 *
## trans.htrans.hal:treatmentdry 0.21094 0.24311 0.868 0.3884
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3844 on 72 degrees of freedom
## Multiple R-squared: 0.8429, Adjusted R-squared: 0.8276
## F-statistic: 55.19 on 7 and 72 DF, p-value: < 2.2e-16
#both trans, no cis
test<-sim.ct(wts=c(0,2,2,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00672 -0.21342 -0.00887 0.27428 0.66463
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.02468 0.20514 24.494 <2e-16 ***
## cisfil 0.03623 0.16750 0.216 0.8294
## trans.ftrans.fil 2.07165 0.16750 12.368 <2e-16 ***
## trans.htrans.hal 1.84003 0.16750 10.985 <2e-16 ***
## treatmentdry -0.21039 0.29011 -0.725 0.4707
## cisfil:treatmentdry 0.14065 0.23688 0.594 0.5545
## trans.ftrans.fil:treatmentdry -0.23421 0.23688 -0.989 0.3261
## trans.htrans.hal:treatmentdry 0.50802 0.23688 2.145 0.0354 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3745 on 72 degrees of freedom
## Multiple R-squared: 0.8918, Adjusted R-squared: 0.8812
## F-statistic: 84.74 on 7 and 72 DF, p-value: < 2.2e-16
#trade off
test<-sim.ct(wts=c(-1,1,-1,0,0,0,0), nind=80, type="normal")
##
## Call:
## lm(formula = value ~ 1 + cis + trans.f + trans.h + cis * treatment +
## trans.f * treatment + trans.h * treatment, data = grps.long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.92377 -0.21170 -0.02124 0.26371 0.90692
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.01068 0.19925 25.148 < 2e-16 ***
## cisfil -1.11749 0.16269 -6.869 1.93e-09 ***
## trans.ftrans.fil 1.10666 0.16269 6.802 2.55e-09 ***
## trans.htrans.hal -1.05469 0.16269 -6.483 9.78e-09 ***
## treatmentdry -0.06833 0.28178 -0.242 0.809
## cisfil:treatmentdry 0.16533 0.23007 0.719 0.475
## trans.ftrans.fil:treatmentdry -0.05253 0.23007 -0.228 0.820
## trans.htrans.hal:treatmentdry 0.08254 0.23007 0.359 0.721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3638 on 72 degrees of freedom
## Multiple R-squared: 0.6982, Adjusted R-squared: 0.6689
## F-statistic: 23.8 on 7 and 72 DF, p-value: < 2.2e-16
