-
Notifications
You must be signed in to change notification settings - Fork 0
/
final negative grouped.R
175 lines (132 loc) · 5.31 KB
/
final negative grouped.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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
#Final negative groups
#load data
care=read.table(file="caredata.csv", header=TRUE, sep=",")
names(care)
str(care)
#Installed Packages
install.packages("lme4")
install.packages("geepack")
install.packages("lsmeans")
install.packages("dplyr")
install.packages("MuMIn")
install.packages("ggplot2")
install.packages("nortest")
install.packages("multcomp")
install.packages("Matrix")
#Library to run
library(nlme)
library(dplyr)
#packages gotten from last years project
library(lattice)
library(lme4)
library(geepack)
library(lsmeans)
library(dplyr)
library(MuMIn)
library(ggplot2)
library(nortest)
library(multcomp)
#Putting the 4 condition in as a factor
care$cond.f=as.factor(care$cond)
str(care)
#Form the negative count
care$negativect=rowSums(care[,c("iproxct", "iproxout","threatg", "bark", "headshake", "yawn", "aggb", "pace", "aggdisp", "teeth", "rubgen", "sway", "turn", "scream")], na.rm=T)
care$negativect
#Negative count behaviors modelling
m1=gls(negativect~cond.f, data=care, na.action=na.omit, method="ML")
summary(m1)
#model2 - try put individual in as a main factor so cond.f+id
m2=lme(negativect~cond.f, random=~1|id,data=care, na.action=na.omit, method="ML")
summary(m2)
anova(m1,m2)
#Since the model with the random factor had a significant lower IAC score we wanna use that.
#M2 AIC = 1731.47 It is significant P = 0,0141
#Plotting residuals for m2 to check if we can use this model
op=par(mfrow=c(2,2), mar=c(5,4,1,2))
plot(m2, add.smooth=FALSE, which=1)
E=resid(m2)
hist(E,xlab="residuals", main="")
plot(filter(care, !is.na(negativect)) %>% dplyr::select(id), #this code will filter NAs out. it's best to use this nested in the plot command all the time
E, xlab="Treatment", ylab="residuals")
plot(filter(care, !is.na(negativect)) %>% dplyr::select(cond.f), #this code will filter NAs out. it's best to use this nested in the plot command all the time
E, xlab="Treatment", ylab="residuals")
qqnorm(residuals(m2))
qqline(residuals(m2))
ad.test(residuals(m2))#This ad.test is also significant P=2.2*10^-16
summary(m2)
#Now we test for differences between condition
odel.matrix.gls <- function(m2, ...){
model.matrix(terms(m2), data = getData(m2), ...)
}
model.frame.gls <- function(m2, ...){
model.frame(formula(m2), data = getData(m2), ...)
}
terms.gls <- function(m2, ...){
terms(model.frame(m2),...)
}
multCompTukey <- glht(m2, linfct = mcp(cond.f = "Tukey"))
summary(multCompTukey)
#There are significante between 2-1, 3-1, 4-1.
lsmeans(m2,pairwise~cond.f)#This one schould be able to show how the behaviors have changes
#There are significant difference between 1-2 P=0,0058, 1-3 P=0.0002, 1-4 P=0.0017.
#We now normalize the data
care$l.negativect=log(care$negativect+1)
#Now we run same models with the normalised data
m1=gls(l.negativect~cond.f, data=care, na.action=na.omit, method="ML")
summary(m1)
#model2 - try put individual in as a main factor so cond.f+id
m2n=lme(l.negativect~cond.f, random=~1|id,data=care, na.action=na.omit, method="ML")
summary(m2n)
anova(m1,m2n)
# Model 2 is still significant different.
op=par(mfrow=c(2,2), mar=c(5,4,1,2))
plot(m2n, add.smooth=FALSE, which=1)
E=resid(m2n)
hist(E,xlab="residuals", main="")
plot(filter(care, !is.na(l.negativect)) %>% dplyr::select(id), #this code will filter NAs out. it's best to use this nested in the plot command all the time
E, xlab="Treatment", ylab="residuals")
plot(filter(care, !is.na(l.negativect)) %>% dplyr::select(cond.f), #this code will filter NAs out. it's best to use this nested in the plot command all the time
E, xlab="Treatment", ylab="residuals")
qqnorm(residuals(m2n))
qqline(residuals(m2n))
ad.test(residuals(m2n))#
summary(m2n)
#Now we test for differences between condition
odel.matrix.gls <- function(m2n, ...){
model.matrix(terms(m2n), data = getData(m2n), ...)
}
model.frame.gls <- function(m2n, ...){
model.frame(formula(m2n), data = getData(m2n), ...)
}
terms.gls <- function(m2n, ...){
terms(model.frame(m2n),...)
}
multCompTukey <- glht(m2n, linfct = mcp(cond.f = "Tukey"))
summary(multCompTukey)
#There are significant results between 1-2, 1-3, 1-4,
lsmeans(m2n,pairwise~cond.f)
#With this one there are 3 significant results.Same as before.
#Test variance equality if the p is sign run alternate variance structure
bartlett.test(l.negativect~id, data=care)
#P is not significant so alternate variance structure is not needed.
#Now we see if alternate variant structures help
#Now we run the alternate variance structures for the normalised model
vf1=varIdent(form = ~1|id)
vf2=varIdent(form = ~1|cond.f)
#Now using these on the normalized model
#alternate models
M2n=lme(l.negativect~cond.f, random=~1|id,data=care, na.action=na.omit, method="ML")
M2.1<-lme(l.negativect~cond.f, random=~1|id,data=care, na.action=na.omit, weights=vf1)
M2.2<-lme(l.negativect~cond.f, random=~1|id,data=care, na.action=na.omit, weights=vf2)
#checking which model is best
anova(M2.2,M2.1)
#Since the normalised model without the alternate structure is better we just finally check for autocorrelation.
#Plotting x with the model
x<-care$l.negativect[!is.na(care$l.negativect)]#removes na values from column
E2<-residuals(M2n,type="normalized")
plot(x, E2)
#residuals v observation
#check for autocorrelation
acf(E2, na.action=na.pass,
main="Auto-correlation plot for residuals")
#There are no autocorrelation.