CODEWITHSUNDEEP
rplotemmeansglmmtmb
Mondik
Mondik

Reputation: 11

In R, plot emmeans of glmmTMB linear model. Error in linkinv(summ$the.emmean) : could not find function "linkinv"

I would like to plot the emmeans of a glmmTMB model using plot(). When my glmmTMB model takes in log transformed data, plot(emmeans(glmmTMB_model)) works just fine. However, when I attempt to plot the emmeans of the glmmTMB model of non-transformed data, I get the following error: Error in linkinv(summ$the.emmean) : could not find function "linkinv".

See my code below:

Site <- c(4,4,5,5)
Treatment <- c("Burnt_Intact-Vegetation","Burnt_Intact-Vegetation", "Burnt_Cleared", "Burnt_Cleared")
pH <- c(5.94, 5.91, 5.44, 5.49)

pH_EC_data <- data_frame(Site, Treatment, pH)

pH_model <- glmmTMB(pH~Treatment+(1|Site), data = pH_EC_data)
log10pH_model<- glmmTMB(log10(`pH`)~Treatment+(1|Site), data = pH_EC_data)


log10pH_analysis <- emmeans(log10pH_model, pairwise~Treatment, type = "response")
plot(log10pH_analysis)
##This plot works just fine.

pH_analysis <- emmeans(pH_model, pairwise~Treatment, type = "response")
plot(pH_analysis)
##This code results in the following error: Error in linkinv(summ$the.emmean) : could not find function "linkinv"

Note, log10pH_analysis and pH_analysis differ by one column. Emmeans of glmmTMB of logged data creates a "response" column whereas the same manipulation of non-tranformed data resulted in an "emmeans" column. See below:

log10pH_analysis

$emmeans
 Treatment               response     SE df lower.CL upper.CL
 Burnt_Cleared               5.47 0.0167  2     5.40     5.54
 Burnt_Intact-Vegetation     5.90 0.0180  2     5.82     5.98

Confidence level used: 0.95 
Intervals are back-transformed from the log10 scale 

$contrasts
 contrast                                  estimate      SE df t.ratio p.value
 Burnt_Cleared - (Burnt_Intact-Vegetation)  -0.0329 0.00188  2 -17.520 0.0032 

Note: contrasts are still on the log10 scale

pH_analysis

$emmeans
 Treatment               emmean     SE df lower.CL upper.CL
 Burnt_Cleared             5.47 0.0176  2     5.39     5.55
 Burnt_Intact-Vegetation   5.90 0.0176  2     5.82     5.98

Confidence level used: 0.95 

$contrasts
 contrast                                  estimate     SE df t.ratio p.value
 Burnt_Cleared - (Burnt_Intact-Vegetation)    -0.43 0.0249  2 -17.238 0.0033

Thank you.

Upvotes: 1

Views: 688

Answers (2)

Russ Lenth
Russ Lenth

Reputation: 6790

I am able to duplicate this with a similar example. It is due to a coding error that I will correct in the next update of the package.

The problem is related to the fact that your pH_analysis object is actually a list of two emmGrid objects -- one for the estimated marginals means, and the other for the pairwise comparisons of them. If you do

plot(pH_analysis, which = 1)    # or plot(pH_analysis[[1]])

it will work just fine. The error occurs in trying to plot the second one as well, and the coding error has to do with the fact that the pairwise comparisons are not on a transformed scale but it thinks they are.

I do suggest that you not use plot() on the results of an emmeans() call with a two-sided formula or list of specs; and instead pick out the part you actually want to plot. Also, I think you must have a fairly old version of emmeans, as the default for plot() in recent versions is which = 1.

Thanks for reporting this bug.

Upvotes: 1

Mondik
Mondik

Reputation: 11

It seems that type = "response" causes the problem for the non-transformed data.

pH_analysis <- emmeans(pH_model, pairwise~Treatment)
plot(pH_analysis)
##Plots beautifully; problem resolved

Upvotes: 0

Related Questions