Dennis Kim
Reputation: 13
Prediction using a ns spline fit (glmer with ns spline)
I am using the glmer() function to determine whether invasive earthworm species weights have significant impacts on 4 different tree species in Minnesota. I am trying to predict the probability of events (because it's logistic regression/ whether tree species will survive or not) at a certain earthworm weight after adjusting other variables. We are also thinking of using odds ratio to compare whether there is a difference among different tree species. So I want to use the predict() function to specify a certain weight value and predict the probability of events.
Here is the summary(data). We are specifically looking at meanwormwt (average worm weight). What I ultimately want is to predict the probability of events given a certain mean worm weight value. For example, "what's the probability of the event when the mean worm weight is 0.3246?"
Exclose DensF Species Park Parkclust ParkPlot
No :2015 Min. : 8.00 ACSA:1011 BRSP:672 BRSP1 : 168 BRSP11 : 96
Yes:2016 1st Qu.:28.00 QUMA:1004 CSP :672 BRSP2 : 168 BRSP2 : 96
Median :48.00 RHCA:1009 GLSP:672 BRSP3 : 168 BRSP6 : 96
Mean :37.52 TIAM:1007 GRF :672 BRSP4 : 168 BRSP8 : 96
3rd Qu.:48.00 NLP :671 CAM1 : 168 CSP10 : 96
Max. :48.00 SSP :672 CAM2 : 168 CSP2 : 96
(Other):3023 (Other):3455
Unique x y PlantId Cluster
BRSP11C: 48 Min. :1.000 Min. :1.000 NLP8C.1.2 : 2 Min. :1.0
BRSP11E: 48 1st Qu.:2.000 1st Qu.:2.000 SSP10E.2.5 : 2 1st Qu.:2.0
BRSP2C : 48 Median :3.000 Median :4.000 BRSP10C.1.1: 1 Median :3.0
BRSP2E : 48 Mean :2.976 Mean :4.142 BRSP10C.1.2: 1 Mean :2.5
BRSP6C : 48 3rd Qu.:4.000 3rd Qu.:6.000 BRSP10C.1.3: 1 3rd Qu.:3.5
BRSP6E : 48 Max. :6.000 Max. :8.000 BRSP10C.1.4: 1 Max. :4.0
(Other):3743 (Other) :4023
Plot LowLevXBlks DeerFPP13park DeerFPP15park DeerAvgFrac
Min. : 1.000 SSP10E:2 : 9 Min. : 1.20 Min. :0.9167 Min. :0.1654
1st Qu.: 4.000 BRSP11C:1: 8 1st Qu.: 1.60 1st Qu.:1.9167 1st Qu.:0.4878
Median : 7.000 BRSP11C:2: 8 Median : 1.90 Median :2.4167 Median :0.5365
Mean : 6.392 BRSP11C:3: 8 Mean : 8.13 Mean :2.3194 Mean :0.5420
3rd Qu.: 9.500 BRSP11C:4: 8 3rd Qu.:12.80 3rd Qu.:3.0000 3rd Qu.:0.6673
Max. :12.000 BRSP11C:5: 8 Max. :20.90 Max. :3.3333 Max. :0.8500
(Other) :3982
DeerFPP13plot DeerFPP15plot Lt15CnpyOpen Lt13CnpyOpen pHH2013
Min. : 0.00 Min. :0.000 Min. : 4.04 Min. : 5.01 Min. :5.48
1st Qu.: 0.00 1st Qu.:1.000 1st Qu.: 9.81 1st Qu.: 8.79 1st Qu.:6.46
Median : 1.00 Median :1.500 Median :13.92 Median :13.73 Median :6.89
Mean : 2.29 Mean :1.663 Mean :15.44 Mean :18.27 Mean :6.93
3rd Qu.: 3.00 3rd Qu.:2.000 3rd Qu.:18.51 3rd Qu.:23.69 3rd Qu.:7.48
Max. :11.00 Max. :7.000 Max. :58.81 Max. :57.38 Max. :8.17
NA's :96
worm16ct worm16wt meanwormct meanwormwt PlotMoist
Min. : 0.0 Min. :0.0000 Min. : 0.00 Min. :0.0000 Min. :0.1206
1st Qu.:10.0 1st Qu.:0.2729 1st Qu.:13.00 1st Qu.:0.3246 1st Qu.:0.1622
Median :22.0 Median :0.6463 Median :23.33 Median :0.6292 Median :0.1873
Mean :26.5 Mean :1.1226 Mean :27.92 Mean :0.9185 Mean :0.1948
3rd Qu.:37.0 3rd Qu.:1.2940 3rd Qu.:40.00 3rd Qu.:1.0901 3rd Qu.:0.2239
Max. :94.0 Max. :7.5096 Max. :85.33 Max. :3.6966 Max. :0.2927
NA's :96 NA's :96
ClustMoist AvgJJA AvgDJF TotalBA FireObsF
Min. :0.1366 Min. :18.90 Min. :-9.300 Min. :0.4269 FireObs : 96
1st Qu.:0.1710 1st Qu.:20.90 1st Qu.:-7.200 1st Qu.:1.1455 NotFireObs:3935
Median :0.1883 Median :21.30 Median :-6.600 Median :1.4635
Mean :0.1966 Mean :21.28 Mean :-6.688 Mean :1.5469
3rd Qu.:0.2190 3rd Qu.:21.80 3rd Qu.:-6.300 3rd Qu.:1.9987
Max. :0.2735 Max. :24.50 Max. :-4.700 Max. :2.7619
Alive12Sum Alive12SumF Alive16JuneF Alive16June
Min. :0.0000 Alive:3927 Alive:1556 Min. :0.0000
1st Qu.:1.0000 Dead : 104 Dead :2379 1st Qu.:0.0000
Median :1.0000 NA's : 96 Median :0.0000
Mean :0.9742 Mean :0.3954
3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :1.0000 Max. :1.0000
NA's :96
I used ns spline on mean earthworm weight and divided into 3 splines. What codes do I need to use for this? I've tried to use predict.ns or predict.merMod but I don't know how to since we are not looking for just overall average predict value.. We want to look at predict value at certain weight. What commands should I try? What should I do?
Here is my glmer code:
```{r}
nsglm<-glmer(Mort16JuneAPF ~ Exclose*Species + ns(meanwormwt, df=3, knots=c(0.3246,1.0901))*Species + (1 | Park) + (1 | Cluster:Park) + (1 | Plot:Cluster:Park) + (1|Exclose:ParkPlot) + (1 | x:Unique), data = mydata, family = binomial, control=glmerControl(optimizer="bobyqa", calc.derivs = FALSE))
summary(nsglm)
```
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: Mort16JuneAPF ~ Exclose * Species + ns(meanwormwt, df = 3, knots = c(0.3246,
1.0901)) * Species + (1 | Park) + (1 | Cluster:Park) + (1 |
Plot:Cluster:Park) + (1 | Exclose:ParkPlot) + (1 | x:Unique)
Data: mydata
Control: glmerControl(optimizer = "bobyqa", calc.derivs = FALSE)
AIC BIC logLik deviance df.resid
4253.3 4410.2 -2101.7 4203.3 3910
Scaled residuals:
Min 1Q Median 3Q Max
-4.2758 -0.6492 0.2821 0.6346 4.0010
Random effects:
Groups Name Variance Std.Dev.
x:Unique (Intercept) 0.01345 0.1160
Exclose:ParkPlot (Intercept) 0.51799 0.7197
Plot:Cluster:Park (Intercept) 0.00000 0.0000
Cluster:Park (Intercept) 0.28753 0.5362
Park (Intercept) 0.03863 0.1965
Number of obs: 3935, groups:
x:Unique, 564; Exclose:ParkPlot, 142; Plot:Cluster:Park, 71; Cluster:Park, 24; Park, 6
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6957 0.4435 1.569 0.116746
ExcloseYes -2.7133 0.2090 -12.981 < 2e-16 ***
SpeciesQUMA 1.2551 0.3827 3.279 0.001041 **
SpeciesRHCA -0.6303 0.3407 -1.850 0.064331 .
SpeciesTIAM -0.5476 0.3500 -1.565 0.117687
ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))1 1.2171 0.6496 1.874 0.060986 .
ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))2 0.8967 0.9645 0.930 0.352534
ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))3 -0.2013 0.7047 -0.286 0.775132
ExcloseYes:SpeciesQUMA 1.5177 0.2375 6.391 1.65e-10 ***
ExcloseYes:SpeciesRHCA 2.2524 0.2138 10.533 < 2e-16 ***
ExcloseYes:SpeciesTIAM 1.0164 0.2295 4.430 9.44e-06 ***
SpeciesQUMA:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))1 -0.3065 0.6130 -0.500 0.617043
SpeciesRHCA:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))1 -1.0661 0.5614 -1.899 0.057555 .
SpeciesTIAM:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))1 0.6600 0.6074 1.087 0.277240
SpeciesQUMA:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))2 -2.1818 0.8225 -2.653 0.007984 **
SpeciesRHCA:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))2 -1.3299 0.7390 -1.800 0.071897 .
SpeciesTIAM:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))2 3.0146 0.7774 3.878 0.000105 ***
SpeciesQUMA:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))3 -2.8120 0.5579 -5.041 4.64e-07 ***
SpeciesRHCA:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))3 -0.4749 0.5100 -0.931 0.351807
SpeciesTIAM:ns(meanwormwt, df = 3, knots = c(0.3246, 1.0901))3 2.4477 0.5762 4.248 2.16e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Upvotes: 1
Views: 811
Answers (1)
Ben Bolker
Reputation: 226192
You need to use predict() with the newdata argument. You need to specify some value for every fixed-effect input variable, e.g.
nd <- with(mydata,
expand.grid(Exclose=levels(Exclose), Species=levels(Species))
nd$meanwormwt <- 0.361
predict(nsglm, re.form=~0, newdata=nd)
re.form=~0 specifies that you want to make population-level predictions (i.e., for a new/unknown value of the random effects grouping factors).
Upvotes: 2
Related Questions
- Using the ns() function from the splines package inside glm()
- predict newdata with NA for smooth.splines
- R - Extract ns spline object from lmer model and predict on new data
- How to predict from ns spline parameters without model object
- fit a smoothing spline using gam function
- Extracting predictions from a GAM model with splines and lagged predictors
- Using spline model in Raster Predict
- Prediction using a natural spline fit
- using ns() for spline prediction
- Only plotting the fitted spline line and not the data points