0
$\begingroup$

I have some data for carbon assimilation vs tree size for a range of tree species. I'm running a mixed model analysis using lme4 in R, with the random effect being species. Ultimately what I'd like to do is simplify the model outputs in such a way that the regression coefficients can be used in an Excel sheet to input the independent variable (tree size) in one cell, and have the dependent variable (carbon assimilation) output in another.

But, the outputs I'm getting from the model for intercept and slope coefficients summary command don't seem to stack up and I wonder if I'm missing something.

Here's a sample of the data along with some code.

library(ggplot2)
library(lme4)
library(merTools)

DBH<-c(4.13, 6.74, 9.35, 11.96, 14.57, 17.18, 19.79, 22.4, 25.01, 27.62, 30.23, 32.84, 35.45, 38.06, 40.67, 43.28, 45.89, 48.5, 51.11, 53.72, 56.33, 58.94, 61.55, 64.16, 66.77, 69.38, 71.99, 74.6, 77.21, 79.82, 82.43, 85.04, 87.65, 90.26, 92.87, 95.48, 98.09, 100.7, 103.31, 105.92, 108.53, 111.14, 113.75, 116.36, 118.97, 121.58, 124.19, 126.8, 129.41, 132.02, 134.63, 137.24, 139.85, 142.46, 145.07, 147.68, 150.29, 152.9, 155.51, 158.12, 160.73, 163.34, 165.95, 168.56, 171.17, 173.78, 176.39, 179, 181.61, 184.22, 186.83, 189.44, 192.05, 194.66, 197.27, 199.88, 202.49, 205.1, 207.71, 210.32, 212.93, 215.54, 218.15, 220.76, 223.37, 225.98, 228.59, 231.2, 233.81, 236.42, 239.03, 241.64, 244.25, 246.86, 249.47, 252.08, 254.69, 257.3, 259.91, 262.52, 265.13, 2.65, 4.04, 5.43, 6.82, 8.21, 9.6, 10.99, 12.38, 13.77, 15.16, 16.55, 17.94, 19.33, 20.72, 22.11, 23.5, 24.89, 26.28, 27.67, 29.06, 30.45, 31.84, 33.23, 34.62, 36.01, 37.4, 38.79, 40.18, 41.57, 42.96, 44.35, 45.74, 47.13, 48.52, 49.91, 51.3, 52.69, 54.08, 55.47, 56.86, 58.25, 59.64, 61.03, 62.42, 63.81, 65.2, 66.59, 67.98, 69.37, 70.76, 72.15, 73.54, 74.93, 76.32, 77.71, 79.1, 80.49, 81.88, 83.27, 84.66, 86.05, 87.44, 88.83, 90.22, 91.61, 93, 94.39, 95.78, 97.17, 98.56, 99.95, 101.34, 102.73, 104.12, 105.51, 106.9, 108.29, 109.68, 111.07, 112.46, 113.85, 115.24, 116.63, 118.02, 119.41, 120.8, 122.19, 123.58, 124.97, 126.36, 127.75, 129.14, 130.53, 131.92, 133.31, 134.7, 136.09, 137.48, 138.87, 140.26, 141.65, 5.65, 8.26, 10.87, 13.48, 16.09, 18.7, 21.31, 23.92, 26.53, 29.14, 31.75, 34.36, 36.97, 39.58, 42.19, 44.8, 47.41, 50.02, 52.63, 55.24, 57.85, 60.46, 63.07, 65.68, 68.29, 70.9, 73.51, 76.12, 78.73, 81.34, 83.95, 86.56, 89.17, 91.78, 94.39, 97, 99.61, 102.22, 104.83, 107.44, 110.05, 112.66, 115.27, 117.88, 120.49, 123.1, 125.71, 128.32, 130.93, 133.54, 136.15, 138.76, 141.37, 143.98, 146.59, 149.2, 151.81, 154.42, 157.03, 159.64, 162.25, 164.86, 167.47, 170.08, 172.69, 175.3, 177.91, 180.52, 183.13, 185.74, 188.35, 190.96, 193.57, 196.18, 198.79, 201.4, 204.01, 206.62, 209.23, 211.84, 214.45, 217.06, 219.67, 222.28, 224.89, 227.5, 230.11, 232.72, 235.33, 237.94, 240.55, 243.16, 245.77, 248.38, 250.99, 253.6, 256.21, 258.82, 261.43, 264.04, 266.65, 4.14, 6.14, 8.14, 10.14, 12.14, 14.14, 16.14, 18.14, 20.14, 22.14, 24.14, 26.14, 28.14, 30.14, 32.14, 34.14, 36.14, 38.14, 40.14, 42.14, 44.14, 46.14, 48.14, 50.14, 52.14, 54.14, 56.14, 58.14, 60.14, 62.14, 64.14, 66.14, 68.14, 70.14, 72.14, 74.14, 76.14, 78.14, 80.14, 82.14, 84.14, 86.14, 88.14, 90.14, 92.14, 94.14, 96.14, 98.14, 100.14, 102.14, 104.14, 106.14, 108.14, 110.14, 112.14, 114.14, 116.14, 118.14, 120.14, 122.14, 124.14, 126.14, 128.14, 130.14, 132.14, 134.14, 136.14, 138.14, 140.14, 142.14, 144.14, 146.14, 148.14, 150.14, 152.14, 154.14, 156.14, 158.14, 160.14, 162.14, 164.14, 166.14, 168.14, 170.14, 172.14, 174.14, 176.14, 178.14, 180.14, 182.14, 184.14, 186.14, 188.14, 190.14, 192.14, 194.14, 196.14, 198.14, 200.14, 202.14, 204.14)
CARBON<-c(0, 0, 10, 20, 30, 50, 70, 100, 130, 160, 210, 250, 310, 370, 440, 510, 590, 680, 770, 880, 990, 1110, 1230, 1370, 1510, 1660, 1820, 1990, 2170, 2350, 2550, 2750, 2970, 3190, 3420, 3670, 3920, 4180, 4460, 4740, 5040, 5340, 5660, 5990, 6330, 6680, 7040, 7410, 7800, 8200, 8600, 9020, 9460, 9900, 10360, 10830, 11310, 11810, 12310, 12830, 13370, 13910, 14470, 15050, 15630, 16230, 16850, 17480, 18120, 18770, 19440, 20130, 20820, 21540, 22260, 23000, 23760, 24530, 25320, 26120, 26930, 27760, 28610, 29470, 30350, 31240, 32150, 33070, 34010, 34960, 35930, 36920, 37920, 38940, 39980, 41030, 42100, 43180, 44290, 45400, 46540, 0, 0, 0, 10, 10, 10, 20, 20, 30, 40, 50, 60, 70, 90, 100, 120, 140, 160, 180, 210, 230, 260, 290, 320, 350, 390, 430, 470, 510, 550, 600, 650, 700, 750, 800, 860, 920, 980, 1050, 1120, 1190, 1260, 1340, 1410, 1490, 1580, 1660, 1750, 1840, 1940, 2040, 2140, 2240, 2350, 2460, 2570, 2680, 2800, 2920, 3050, 3180, 3310, 3440, 3580, 3720, 3860, 4010, 4160, 4310, 4470, 4630, 4800, 4960, 5140, 5310, 5490, 5670, 5860, 6040, 6240, 6430, 6630, 6840, 7040, 7260, 7470, 7690, 7910, 8140, 8370, 8600, 8840, 9080, 9330, 9580, 9830, 10090, 10350, 10620, 10890, 11170, 0, 0, 10, 10, 20, 30, 40, 50, 70, 80, 100, 120, 140, 160, 180, 210, 230, 260, 290, 310, 340, 370, 400, 420, 450, 490, 520, 550, 580, 620, 650, 690, 730, 770, 800, 840, 880, 920, 970, 1010, 1050, 1100, 1140, 1190, 1240, 1280, 1330, 1380, 1430, 1480, 1530, 1590, 1640, 1690, 1750, 1800, 1860, 1910, 1970, 2030, 2090, 2150, 2210, 2270, 2330, 2400, 2460, 2520, 2590, 2650, 2720, 2790, 2860, 2920, 2990, 3060, 3130, 3200, 3280, 3350, 3420, 3500, 3570, 3650, 3720, 3800, 3880, 3960, 4040, 4120, 4200, 4280, 4360, 4440, 4530, 4610, 4700, 4780, 4870, 4950, 5040, 0, 0, 10, 20, 30, 40, 50, 70, 90, 110, 130, 160, 190, 230, 270, 310, 350, 400, 460, 510, 570, 640, 700, 780, 850, 930, 1020, 1110, 1200, 1300, 1400, 1510, 1620, 1740, 1860, 1990, 2120, 2250, 2390, 2540, 2690, 2850, 3010, 3170, 3340, 3520, 3700, 3890, 4080, 4280, 4480, 4690, 4910, 5130, 5350, 5590, 5820, 6070, 6320, 6570, 6830, 7100, 7370, 7650, 7940, 8230, 8520, 8830, 9140, 9450, 9770, 10100, 10440, 10780, 11120, 11480, 11840, 12200, 12580, 12960, 13340, 13740, 14140, 14540, 14960, 15380, 15800, 16240, 16680, 17120, 17580, 18040, 18510, 18980, 19460, 19950, 20450, 20950, 21460, 21980, 22500)
SPECIES<-rep(c("OAK", "GINKGO","SWEETGUM", "MAGNOLIA"), each = 101)

data<-data.frame(DBH,CARBON,SPECIES)
str(data)

### PLOT THE DATA BY SPECIES

ggplot(data,aes(y=CARBON, x=DBH,col=SPECIES))+geom_point()+
  geom_line()+facet_wrap((~SPECIES))

## RUN THE MODEL WITH A SECOND ORDER POLYNOMIAL 

model<-lmer(CARBON ~ 1 + poly(DBH,2) +(1|SPECIES)+(0+poly(DBH,2)|SPECIES),data=data)

## IF IT FAILS TO CONVERGE, RUN IT AGAIN WITH THE PARAMETER ESTIMATES PRIOR TO FAILED CONVERGENCE

    ss <- getME(model,c("theta","fixef"))
    model <- update(model,start=ss)

summary(model)

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)   7.299e+03  1.993e+03 3.001e+00   3.663   0.0352 *
poly(DBH, 2)1 1.622e+05  4.588e+04 3.003e+00   3.536   0.0384 *
poly(DBH, 2)2 5.804e+04  1.787e+04 3.005e+00   3.248   0.0474 *

# GENERATE PREDICTED VALUES FOR EACH RANDOM EFFECT

newdat<-expand.grid(DBH=seq(0,300,length=602),
                  SPECIES=levels(data$SPECIES))
newdat$CARBON<-predict(model,newdata=newdat)

# GENERATE PREDICTED VALUES FOR THE MEAN

newdat2<-data.frame(DBH=seq(0,300,length=602))
newdat2$CARBON<-predict(model,newdata=newdat2,re.form=~0)

# PLOT THE DATA


plot<-ggplot(data,aes(x=DBH,y=CARBON,group=SPECIES,col=SPECIES))+
  geom_line(data=newdat,lwd=0.3)+
   geom_line(data=newdat2,aes(group=NULL),col="black",lwd=1)
plot

## ADD IN THE PREDICTED VALUES FOR THE MEAN AS POINTS, TO CONFIRM THEY'RE IN THE RIGHT PLACE

plot+geom_point(data=newdat2,aes(group=NULL),col="red",cex=0.4)


# RUN A SIMPLE POLYNOMIAL REGRESSION ON THE newdat2 DATA - THE DATA POINTS FOR THE MEAN FIT LINE


lm_mod<-lm(CARBON ~ poly(DBH,2), newdat2)
    summary(lm_mod)

Coefficients:
               Estimate Std. Error   t value Pr(>|t|)    
(Intercept)   1.357e+04  5.770e-13 2.351e+16   <2e-16 ***
poly(DBH, 2)1 3.091e+05  1.416e-11 2.184e+16   <2e-16 ***
poly(DBH, 2)2 9.652e+04  1.416e-11 6.819e+15   <2e-16 ***

The slope and intercept coefficients don't seem to fit the shape of the mean regression curve. So, I export the data to Excel

setwd("C:/YOUR FILE NAME")
write.csv(newdat2, file="DUMMY_DATA.csv", row.names = F)

Now when I insert a scatter plot and add a second order polynomial trendline, I get the following:

y = 0.5845x^2 - 30.104x + 531.98

The slope and intercept coefficients for this curve match the data which R exported, so why do they differ from the slope and intercept coefficients in R?

I can then use these coefficients in a couple of cells to generate an easy formula to check the correctness, which seems to stack up.

The overall objective of the analysis is to have the outputs from R simplified so that an unfamiliar person can use Excel to simply input a few parameters to get an output based on the analyses, using the mean regression equation from the mixed model. But, why don't the coefficients I'm getting in R seem right, and why don't they match those from Excel when I use the same data? What am I missing?

$\endgroup$
  • 1
    $\begingroup$ Which coefficients don't seem right? Can you revise your post by tying explicitly the coefficients to the models they come from? For example: "the coefficients from the lm model don't seem right". $\endgroup$ – Isabella Ghement Feb 16 at 16:24
  • 1
    $\begingroup$ Hi @IsabellaGhement. I've updated the question. Specifically, what I need to know is why don't the slope and intercept coefficients in R for the mean regression curve from the mixed model match the slope and intercept coefficients which Excel gives me when I generate a set of fitted values from the mixed model? The ones in Excel seem correct and to fit the data. $\endgroup$ – A.Benson Feb 18 at 22:25
  • $\begingroup$ Thanks! See my answer. $\endgroup$ – Isabella Ghement Feb 20 at 14:06
  • $\begingroup$ Thanks for accepting my answer. For completeness, I edited my original answer after you accepted it to include some R output and an Addendum. $\endgroup$ – Isabella Ghement Feb 20 at 20:18
  • $\begingroup$ If you fit the model<-lmer(CARBON ~ 1 + poly(DBH,2) +(1|SPECIES)+(0+poly(DBH,2)|SPECIES),data=data), you can add the option raw = TRUE inside all your poly() terms to force the model to fit a raw second-degree polynomial of order 2 for each tree species included in the data. Once you fit the model - presuming it does not give you a singular fit - you can just use fixef(model) to extract the fixed effects associated with (raw) DBH and squared (raw) DBH for the "typical" species. Is this what you are looking for? $\endgroup$ – Isabella Ghement Mar 10 at 3:01
2
$\begingroup$

@A.Benson: Thank you for the updates to your post. When you state "Now when I insert a scatter plot and add a second order polynomial trendline, I get the following:", it sounds like you do this in Excel. Specifically, you use Excel to fit a quadratic polynomial fit to the scatterplot of CARBON versus DBH using the data newdat2.

While I don't use Excel for model fitting, I suspect that Excel uses the raw powers of DBH for the model, instead of orthogonal polynomial terms. In other words, Excel uses the equivalent of the following in R:

lm_mod_raw <- lm(CARBON ~ poly(DBH,2, raw = TRUE), newdat2)

coef(lm_mod_raw)

which is the same as:

 lm_mod_raw <- lm(CARBON ~ DBH + I(DBH^2), newdat2)

 coef(lm_mod_raw)

However, you are comparing the raw polynomial model fit coefficients produced by Excel against the orthogonal polynomial model fit coefficients produced by R via the commands:

lm_mod <- lm(CARBON ~ poly(DBH,2, raw = FALSE), newdat2)

coef(lm_mod)

which are also equivalent with

lm_mod <- lm(CARBON ~ poly(DBH,2), newdat2)

coef(lm_mod)

The coefficients produced by the two model fits are naturally different since the models used different predictors. Specifically, the coefficients for the model lm_mod_raw are:

> coef(lm_mod_raw)
              (Intercept) poly(DBH, 2, raw = TRUE)1 poly(DBH, 2, raw = TRUE)2 
              531.9820795               -30.1038439                 0.5845085 

while the coefficients for the model lm_mod are:

> coef(lm_mod)
               (Intercept) poly(DBH, 2, raw = FALSE)1 poly(DBH, 2, raw = FALSE)2 
                  13566.25                  309145.28                   96524.24 

You can extract the predictors used by R for the two model fits with these commands:

model.matrix(lm_mod_raw)

model.matrix(lm_mod)

This will enable you to see that the two models use different predictors.

Addendum:

You can calculate the correlation between the predictors included in each model with these commands:

round(cor(model.matrix(lm_mod_raw)[,2:3]),3)

round(cor(model.matrix(lm_mod)[,2:3]),3)

This way, you can see that the raw powers of DBH used in model lm_mod_raw are highly correlated (corr = 0.968):

> round(cor(model.matrix(lm_mod_raw)[,2:3]),3)
                          poly(DBH, 2, raw = TRUE)1 poly(DBH, 2, raw = TRUE)2
poly(DBH, 2, raw = TRUE)1                     1.000                     0.968
poly(DBH, 2, raw = TRUE)2                     0.968                     1.000

whereas the orthogonal polynomials of DBH used in the model lm_mod are no longer correlated (corr = 0):

> round(cor(model.matrix(lm_mod)[,2:3]),3)
                               poly(DBH, 2, raw = FALSE)1 poly(DBH, 2, raw = FALSE)2
    poly(DBH, 2, raw = FALSE)1                          1                          0
    poly(DBH, 2, raw = FALSE)2                          0                          1 

The vif() function in the car package will detect the collinearity present in the model lm_mod_raw due to the high correlation among the model predictors:

library(car)

vif(lm_mod_raw)

by posting the following output and warning:

> vif(lm_mod_raw)
poly(DBH, 2, raw = TRUE)1 poly(DBH, 2, raw = TRUE)2 
                 15.95037                  15.95037 
Warning message:
In summary.lm(object) : essentially perfect fit: summary may be unreliable

The variance inflation factor (VIF) values are high (larger than the threshold of 10 used by some people) for both predictors.

Contrast this with the VIF values for the model lm_mod:

vif(lm_mod)

which are shown below:

> vif(lm_mod)
poly(DBH, 2, raw = FALSE)1 poly(DBH, 2, raw = FALSE)2 
                     1                          1 
Warning message:
In summary.lm(object) : essentially perfect fit: summary may be unreliable

The "essentially perfect fit" warning still remains but the VIF values are now much smaller. You probably wouldn't want to fit a polynomial regression model to predicted values produced by a mixed effects model (especially one where convergence didn't occur) anyway.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks @Isabella Ghement. Can I use raw=T to get the coefficients from the mixed effects model in the same way? What I'm interested in are the raw=T coefficients of the mean fixed effect curve. How can I obtain these without fitting a polynomial to the predicted values? Thanks for your help $\endgroup$ – A.Benson Mar 7 at 7:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.