Interpretation of glm and HSD.test output

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  shapiro.test(residuals(model.best)) tells me, that the residuals are normally distributed, and here's the plot: 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)


returns this:

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 :)

• you ask quite a few questions here, not all very well connected -- this makes the question somehwat unfocused/broad. Please try to make your questions focus on a specific issue, and then if you have other questions that are not intimately connected to the first, see if you can ask them separately. It's possible to link to earlier questions if you need them for context. Since you already have an answer, I'm not going to close this, but generally it's better to keep one main issue per question. Commented Jul 30, 2016 at 3:00
– oli
Commented Jul 30, 2016 at 15:08

You could try some of these approaches:

1. A GLM consists of 3 components - the conditional distribution of the response variable given values of the independent variables, the linear predictor function, and the link function translating the linear predictor to the expected response variable. You've used the default family - Gaussian with default identity link function, this is the same as a multiple linear regression. Check if this is suitable for your target variable mass.loss.a.pct . Evaluate if other alternative family/link functions would be better suited to this dataset.
2. You haven't mentioned the fitted model's null deviance as compared to the fit deviance, that's a good indicator to examine.
3. The leverage plot looks somewhat worrying, it appears you have too many high leverage points which can result in a poor fit.
4. Try using a validation set to see how well the model performs on held-out samples. Since this is a regression model, evaluate the root mean squared error (RMSE), which is the square root of the mean squared error.
5. Since this dataset shows variables with interaction are significant to explain the response variable, you may be able to use regression trees. Try fitting an rpart model on this dataset.

3. Since these are nominal variables, GLM expands each of these into multiple dummy variables. So, the number of predictors grows to as many as the number of levels of the variable. With interaction effects, these multiply to an even greater number of variables. Your final model is mass.loss.a.pct ~ origin + species + origin:species + 1 . So the number of actual variables used by GLM = Number of unique levels of origin + number of unique levels of species + (number of levels of origins x number of levels of species). Think of each of these dummy variables along with non-dummy variables as a dimension, so if there are 118 dimensions, you're trying to fit a surface through this space. If for any of these dimensions you don't have sufficient values to fit this surface, these coefficients will show up as NA.