Dear Steven,
I am so happy, that you answered me! I tried what you said to put all the
variable in one dataframe. The Pretreatment is not really necessary, because it
didn't show any significance. I didn't copy it into the help, because I tried
to concentrate on the essential things.
Here is the whole code:
PAMdata <-read.table("PAMdata.csv",sep = ";",header=TRUE) #warning: should
have no extra anything in your column names
##inspect##
head(PAMdata)
summary(PAMdata)
str(PAMdata)
##end inspect##
install.packages(c("sciplot","nlme","multcomp"))
library(mvtnorm)
library(splines)
library(survival)
library(sciplot)
library(nlme)
library(multcomp)
#Factors
PAMdata$provenance[PAMdata$provenance == "1"] = "BG"
PAMdata$provenance[PAMdata$provenance == "2"] = "DE"
PAMdata$provenance[PAMdata$provenance == "3"] = "IT"
PAMdata$provenance[PAMdata$provenance == "4"] = "SE"
PAMdata$provenance[PAMdata$provenance == "5"] = "ES"
PAMdata$provenance[PAMdata$provenance == "6"] = "HU"
PAMdata$Treatmentf <- factor(PAMdata$treatment, levels=c("C","F"))
PAMdata$Datef <- factor(PAMdata$Date, levels=c( "25.05.10 14:00","26.05.10
19:00","27.05.2010 7:30","27.05.10 14:00","01.06.10 14:00","02.06.10
19:00","23.06.10 12:30"),ordered=TRUE)
PAMdata$Pretreatmentf <- as.factor(PAMdata$pretreatment)
PAMdata$Provenancef <- as.factor(PAMdata$provenance)
PAMdata$Greenhousef <- as.factor(PAMdata$greenhouse)
PAMdata$Individualf <- as.factor(PAMdata$individual)
PAMdata$PAMval <- (PAMdata$DataPAM)
PAMdata$Code<-(PAMdata$code)
head(PAMdata)
####### Statistischer Test mit Anova #########
summary(PAMaov<-aov(PAMval~Treatmentf*Pretreatmentf*Provenancef+Error(Datef/Code),data
= PAMdata))
#############################################################
##Linear fixed effects model lme
summary(PAM.lme<-lme(PAMval~Treatmentf*Provenancef*Pretreatmentf,
random=~1|Datef/Code, data = PAMdata, na.action=na.omit))
### Tukey test ##
summary(glht(PAM.lme, linfct = mcp(Provenancef = "Tukey")))
Error message:
Fehler in glht.matrix(model = list(modelStruct = list(reStruct = list(Code =
0.808654423456211, :
ncol(linfct) is not equal to length(coef(model))
Zusätzlich: Warnmeldung:
In mcp2matrix(model, linfct = linfct) :
covariate interactions found -- default contrast might be inappropriate
summary(glht(PAM.lme, linfct = mcp(Treatmentf = "Tukey")))
--> gives the same error
traceback()
--> The whole traceback thing is huge, do you really wanna have it? Here are
the last lines:
4: do.call("glht", args)
3: glht.mcp(PAM.lme, linfct = mcp(Provenancef = "Tukey"))
2: glht(PAM.lme, linfct = mcp(Provenancef = "Tukey"))
1: summary(glht(PAM.lme, linfct = mcp(Provenancef = "Tukey")))
with(PAMdata, table(Pretreatmentf, Provenancef, Treatmentf))
--> gives:
, , Treatmentf = C
Provenancef
Pretreatmentf BG DE ES HU IT SE
0 0 0 0 0 0
C 63 63 42 63 63 63
W 63 63 63 63 63 63
, , Treatmentf = F
Provenancef
Pretreatmentf BG DE ES HU IT SE
0 0 0 0 0 0
C 63 63 63 63 63 63
W 63 63 63 63 63 63
sessionInfo()
> sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)
locale:
[1] de_DE.UTF-8/de_DE.UTF-8/C/C/de_DE.UTF-8/de_DE.UTF-8
attached base packages:
[1] splines stats graphics grDevices utils datasets methods base
other attached packages:
[1] multcomp_1.2-4 nlme_3.1-97 sciplot_1.0-7 survival_2.35-8
mvtnorm_0.9-92
loaded via a namespace (and not attached):
[1] grid_2.12.0 lattice_0.19-13
I would like very much to share the data, I am just not exactly sure how to
make it?
Thank you so much for your answer and your help! Now I don't feel so lonely
with my problem anymore!
Best wishes,
Lilith
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
