-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathBGGN273_week9 and 10.R
100 lines (63 loc) · 2.47 KB
/
BGGN273_week9 and 10.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
####Week 9 & 10
diet<-read.csv(file = "DietAncovaData.csv")
#Does type of diet modify the ammount of food consumed?
hist(diet$DietConsumed)
ks.test(diet$DietConsumed,pnorm,mean(diet$DietConsumed),sd(diet$DietConsumed))
model<-lm(DietConsumed~Diet,data=diet)
mean(diet$DietConsumed[diet$Diet=="A"])
mean(diet$DietConsumed[diet$Diet=="B"])
summary(model)
t.test(diet$DietConsumed~diet$Diet)
#Does diet consumed change with increasing mass?
model2<-lm(diet$DietConsumed~diet$Weight)
summary(model2)
#Interaction Model
intmod<-lm(diet$DietConsumed~diet$Diet*diet$Weight)
step(intmod)
summary(intmod)
Bigintmod<-lm(diet$DietConsumed~diet$Diet*diet$Weight*diet$Temperature)
Bigintmod
betterbigmod<-lm(diet$DietConsumed~diet$Diet:diet$Weight+diet$Diet:diet$Temperature+diet$Diet+diet$Weight+diet$Temperature)
step(betterbigmod)
####Randome effects model
library(car)
load("splityield.R")
yields
lm(yield~irrigation,data=yields)
mean(yields$yield)
mean(yields$yield[yields$irrigation=="control"])
mean(yields$yield[yields$irrigation=="irrigated"])
mean(scale(yields$yield,scale=FALSE))
zmodel<-lm(scale(yields$yield,scale=FALSE)~yields$irrigation)
mean(yields$yield)-mean(yields$yield[yields$irrigation=="control"])
str(yields)
comb.field.irrigation<-paste(yields[,2],yields[,3],sep="")
boxplot(yields[,1]~comb.field.irrigation)
library(lme4)
model<-lmer(yield~irrigation+(1|field),data=yields)
nullmodel<-lmer(yield~(1|field),data=yields)
summary(model)
AIC(model)
AIC(nullmodel)
#Gotta have that p value
anova(nullmodel,model)
coef(model) #The intercepts are different but not the slopes! is that okay?
#Random slope
model2<-lmer(yield~irrigation+(1+irrigation|field),data=yields)
coef(model2)
anova(nullmodel,model2)
MegaMODEL<-lmer(yield ~ irrigation * density * fertilizer+ (fertilizer | field) + (density | field)+(irrigation | field),data=yields)
#Model selection in 2 parts
#part 1 decide on your error structure
summary(MegaMODEL)
Errormodel1<-lmer(yield~(fertilizer | field) + (density | field)+(irrigation | field),data=yields)
errorinterceptfull<-lmer(yield~(1|fertilizer) + (1|density)+(1|irrigation ),data=yields)
AIC(errorinterceptfull)
errorinterceptFERT<-lmer(yield~ (1|density)+(1|irrigation ),data=yields)
AIC(errorinterceptFERT)
errorinterceptDIV<-lmer(yield~(1|fertilizer) +(1|irrigation ),data=yields)
AIC(errorinterceptDIV)
errorinterceptIRR<-lmer(yield~(1|fertilizer) + (1|density),data=yields)
AIC(errorinterceptIRR)
JustIRR<-lmer(yield~(1|irrigation ),data=yields)
AIC(JustIRR)