I've posted this on stack overflow a few days ago, but since nobody replied, I assume that it was the wrong forum for this kind of issue. And probably it's way to basic... but as I seriously need some help with this, I re-post it here (with some extensions) and I appreciate any comments or answers.
My datset looks like this:
'data.frame': 124 obs. of 28 variables: $ loglabel : Factor w/ 121 levels "FFBET1B","FFBET2B",..: 62 63 64 65 66 92 93 94 95 95 ... $ origin : Factor w/ 2 levels "F","S": 2 2 2 2 2 2 2 2 2 2 ... $ incubation : Factor w/ 2 levels "F","S": 1 1 1 1 1 2 2 2 2 2 ... $ species : Factor w/ 10 levels "AGR","BET","FEX",..: 1 1 1 1 1 1 1 1 1 1 ... $ plot : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 4 ... $ log : Factor w/ 6 levels "A","B","C","CX",..: 2 1 3 6 5 6 5 2 5 5 ... $ mass.loss.a.pct : num -3.4 5.9 32.5 -0.8 9.4 21.8 11.8 4.8 12.6 NA ...
I am trying to model mass loss from all 5 factors and all possible combinations of these factors.
Firstly, I made an index of rows, which don't contain null:
idx.notNull <- which(rowSums(is.na(Loglife[,c('mass.loss.a.pct','species','incubation', 'origin','plot','log')])) == 0)
My command for the model is:
Lmod1.full <- glm(mass.loss.a.pct ~ species + origin + incubation + plot + log + species:origin + species:incubation + species:plot + species:log + origin:incubation + origin:plot + origin:log + incubation:plot + incubation:log + plot:log, data = Loglife[idx.notNull,],na.action = "na.fail")
After dredging (
model.list <- dredge(Lmod1.full,m.max = 10)) and picking the best model based on df and AICc, the output from
summary(model.best) looks like this:
Call: glm(formula = mass.loss.a.pct ~ origin + species + origin:species + 1, data = Loglife[idx.notNull, ], na.action = "na.fail") Deviance Residuals: Min 1Q Median 3Q Max -41.860 -6.684 -0.950 8.601 29.389 Coefficients: (8 not defined because of singularities) Estimate Std. Error t value Pr(>|t|) (Intercept) 37.602 7.324 5.134 1.42e-06 *** originS -27.091 5.928 -4.570 1.41e-05 *** speciesBET -3.269 8.493 -0.385 0.701138 speciesFEX -14.552 8.383 -1.736 0.085699 . speciesFSY -1.417 8.798 -0.161 0.872425 speciesLKA -5.633 6.082 -0.926 0.356568 speciesPAB -31.682 8.383 -3.779 0.000269 *** speciesPME -8.181 5.928 -1.380 0.170661 speciesPOP -8.312 8.383 -0.992 0.323847 speciesPTR 31.049 5.928 5.238 9.16e-07 *** speciesQRO -6.591 5.928 -1.112 0.268880 originS:speciesBET NA NA NA NA originS:speciesFEX NA NA NA NA originS:speciesFSY NA NA NA NA originS:speciesLKA NA NA NA NA originS:speciesPAB 45.759 8.520 5.371 5.20e-07 *** originS:speciesPME NA NA NA NA originS:speciesPOP NA NA NA NA originS:speciesPTR NA NA NA NA originS:speciesQRO NA NA NA NA
So I don't understand the output from
summary() of the glm.
But first question: Am I right to use glm with this dataset? As you already noticed, I am quite innocent with statistics (but still I'd like to interpret that data I've collected... haha). I started the analysis with mixed models, and due to strange results and someones hint I switched to glm. What I found on the Internet until now, were questions from people with binomial, logistic and Poisson data concerning glm. I can't figure out whether it is ok or not what I've been doing until now.
Assuming glm is fine...
Next question: What does the intercept estimate represent here - 37% mass loss for wood from species AGR from origin F? (the mean of mass loss for AGR at F is slightly different from the intercept estimate; I assume that the difference would be caused by the interaction between species and origin). The next coefficient - does it say that wood (from AGR?) from origin S has 27% lower mass loss than AGR from origin F? For me this would be quite meaningless, as there is no AGR at F. So I wouldn't like to compare every other species to AGR from F, also including the influence of interaction. Is there a possibility to change the intercept in order to make R return a value that has a meaning?
One more: what can I do to get values instead of NAs for the origin:species estimates in the coefficient table? here I found that a possible reason could be that there are more predictors than observations. In my case there are 3 predictors (origin, species, and their interaction) and 111 observations (obtained from
print(idx.notNull)), right? Or is it observation per species? That would be 18 for PAB, and for most others 8 or 9 (also 18 for QRO). I don't see how this is logic?
And does the lack of information about the interaction term reduce the informative value of the model?
Last one: I was told to use
HSD.test() in order to compare the variables.
Loglife.interaction <- cbind(Loglife, origin.spp = interaction(Loglife$origin, Loglife$species)) HSD.test(y = glm(mass.loss.a.pct ~ origin + species + origin.spp, data = Loglife.interaction),'origin.spp',console = T)
Groups, Treatments and means a S.PTR 41.56 a F.FSY 36.19 a F.BET 34.33 a F.QRO 31.01 ab F.POP 29.29 abc S.PAB 24.59 abcd F.FEX 23.05 bcde S.AGR 10.51 cde F.PAB 5.92 cde S.LKA 4.878 de S.QRO 3.92 e S.PME 2.33
I would assume, that R is using the formula from the best model to decide how these groups are split up, so is it already taking the effect of origin and the interaction into consideration when I am looking at only species. But the fact, that origin is listed here makes me insecure about that. Again, I don't what exactly these values mean. So just as an example to make clear what I think they mean: Origin site S and the species specific traits from PTR have different effects on the mass loss in wood. Significantly different to these are the effects from origin site S and the species specific traits from PME on the mass loss in wood (41.56 <-> 2.33). Is that true?
I am convinced I could find all the answers in books and on websites, but actually I am so tired from searching and stumbling upon all those terms which raise even more questions, and also I am running a bit out of time and can't invest much more.
Thanks a lot for reading all this, and I hope there's somebody willing to help me out :)