Reputation: 5
I'm writing this post because I'm stuck in the analysis of a data file from a laboratorial experiment.
In this experiment, I counted the number of females (of a small arthropod) present in a specific environment, throughout 26 time points (TP). However, I want to understand if the number of females was different between each successive time point (e.g. if the number of females counted in TP 1 is different than TP 2; the number of females counted in TP 2 is different than TP 3; and so on...)
The data frame has the following columns:
Replicate (who contain the number of the replicate, going from 1 to 8); TimePoint (the day in which the females where counted, going from 1 to 26); Females (the number of females counted in each time point); and Block (experiment had 2 blocks)
I've tried to do some successive contrasts, but I dont think its the best way. This is my code:
m1<-lmer(Females~TimePoint+(1|Block))
Suc_contrasts2<-glht(m1,linfct=mcp(TimePoint=
c(
"t1 - t2 == 0",
"t2 - t3 == 0",
"t3 - t4 == 0",
"t4 - t5 == 0",
"t5 - t6 == 0",
"t6 - t7 == 0",
"t7 - t8 == 0",
"t8 - t9 == 0",
"t9 - t10 == 0",
"t10 - t11== 0",
"t11 - t12 == 0",
"t12 - t13 == 0",
"t13 - t14 == 0",
"t14 - t15 == 0",
"t15 - t16 == 0",
"t16 - t17 == 0",
"t17 - t18 == 0",
"t18 - t19 == 0",
"t19 - t20 == 0",
"t20 - t21== 0",
"t21 - t22 == 0",
"t22 - t23 == 0",
"t23 - t24 == 0",
"t24 - t25 == 0",
"t25 - t26 == 0")))
summary(Suc_contrasts2)
summary(Suc_contrasts2, test=adjusted ("bonferroni"))
I've been looking on google for other ways to do planned comparisons, but all i've found was not really appropriate for my data set. I'm still new at this, so sorry for any newbie mistake. Thus my question is, is there any better way to compare the number of females I found between each pair of successive time points?
Edit 1:
I also tried doing contrasts like this, but the results don't seem right..
levels(TimePoint)
# [1] "t1" "t10" "t11" "t12" "t13" "t14" "t15" "t16" "t17" "t18" "t19" "t2" "t20" "t21" "t22" "t23" "t24" "t25" "t26"
# [20] "t3" "t4" "t5" "t6" "t7" "t8" "t9"
# tell R which TimePoints to compare
c1<- c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #1v2
c2<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0) #2v3
c3<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0) #3v4
c4<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0) #4v5
c5<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0) #5v6
c5<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0) #6v7
c6<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0) #7v8
c7<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1) #8v9
c8<- c(0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1) #9v10
c9<- c(0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #10v11
c10<- c(0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #11v12
c11<- c(0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #11v12
c12<- c(0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #12v13
c13<- c(0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #13v14
c14<- c(0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #14v15
c15<- c(0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #15v16
c16<- c(0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #16v17
c17<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #17v18
c18<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #18v19
c19<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #19v20
c20<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #20v21
c21<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #21v22
c22<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) #22v23
c23<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0) #23v24
c24<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0) #24v25
c25<- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0, 0, 0) #25v26
# combined the above lines into a matrix
mat <- cbind(c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,c21,c22,c23,c24,c25)
# tell R that the matrix gives the contrasts you want
contrasts(TimePoint) <- mat
model2 <- aov(Females ~ TimePoint)
summary(model2)
# Df Sum Sq Mean Sq F value Pr(>F)
# line2$TimePoint 25 9694 387.8 6.939 <2e-16 ***
# Residuals 390 21794 55.9
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary.aov(model2, split=list(TimePoint=list("1v2"=1, "2v3" = 2, "3v4"=3, "4v5"=4, "5v6"=5, "6v7"=6, "7v8"=7, "8v9"=8, "9v10"=9, "10v11"=10, "11v12"=11, "12v13"=12, "13v14"=13, "14v15"=14, "15v16"=15, "16v17"=16, "17v18"=17, "18v19"=18, "19v20"=19, "20v21"=20, "21v22"=21, "22v23"=22, "23v24"=23, "24v25"=24, "25v26"=25)))
Thanks for your time, André
Upvotes: 0
Views: 369
Reputation: 226057
Another option for fitting successive-differences contrasts:
m1 <- lmer(Females~TimePoint+(1|Block), contrasts=list(TimePoint=MASS::contr.sdif))
This doesn't take the multiplicity of testing into account (which you might get away with since these are pre-planned contrasts): you could use p.adjust()
on the p-values.
@AndreasM's points about the ordering of your factor, choice of random vs fixed effects, etc., should definitely be taken into consideration.
Upvotes: 1
Reputation: 942
I think this website my help you: backward difference coding
Following the information there, difference contrasts between subsequent factor levels could be set like this (see below). Note, that I only use a simple example with 5 factor levels.
#create dummy data
df <- expand.grid(TimePoint = c("t01", "t02", "t03", "t04", "t05"),
Replicate = 1:8, Block = 1:2)
set.seed(2)
df$Females <- runif(nrow(df), min = 0, max = 100)
#set backward difference contrasts
contrasts(df$TimePoint) <- matrix(c(-4/5, 1/5, 1/5, 1/5, 1/5,
-3/5, -3/5, 2/5, 2/5, 2/5,
-2/5, -2/5, -2/5, 3/5, 3/5,
-1/5, -1/5, -1/5, -1/5, 4/5), ncol = 4)
When fitting a simple linear model, the parameter estimates correspond to the expected values, i.e., contrast "TimePoint1" corresponds to t2 - t1, contrast "TimePoint2" to t3 - t2 and so on.
#fit linear model
m1 <- lm(Females ~ TimePoint, data = df)
coef(m1)
(Intercept) TimePoint1 TimePoint2 TimePoint3 TimePoint4
50.295659 -10.116045 7.979465 -10.182389 2.209413
#mean by time point
with(df, tapply(Females, TimePoint, mean))
t01 t02 t03 t04 t05
57.23189 47.11584 55.09531 44.91292 47.12233
I want to add that I am not sure whether what you are trying to do is sensible. But this I don´t feel comfortable to evaluate and it would be a question for CrossValidated. I worry that treating 26 time points as categorical factor levels is not the best way to go. Also, in your initial code you seem to fit a model treating block as a random factor. This does not make sense if block has only 2 levels (as you write), see for example here: Link
Finally, I noticed that in your example, factor levels of your TimePoint variable are not ordered right (t1, t10, t11... instead of t1, t2, t3, ...). You could change this for instance with this line of code:
df$TimePoint <- factor(df$TimePoint, levels = paste0("t", 1:26),
labels = paste0("t", sprintf("%02d", 1:26)))
Upvotes: 0