set.seed(2514)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.0 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.1 ✔ tibble 3.1.8
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(broom.mixed)
## Warning: package 'broom.mixed' was built under R version 4.2.3
library(lme4)
## Warning: package 'lme4' was built under R version 4.2.3
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
data("CO2")
uptake2 = CO2$uptake + runif(84, 0.02, 0.04)
CO2$uptake2 = uptake2
ggplot(CO2, aes(x = conc, y = uptake2, color = Type)) +
geom_point(aes(shape = Treatment)) +
geom_path(aes(group = Plant, lty = Treatment)) +
theme_bw()
Fit1 <- lm(uptake2 ~ I(log(conc)) + Type:Treatment, data = CO2)
Fit2 <- lmer(uptake2 ~ I(log(conc)) + Type:Treatment + (1 | Plant), data = CO2)
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
broom::glance(Fit1)
## # A tibble: 1 × 12
## r.squared adj.r.squa…¹ sigma stati…² p.value df logLik AIC BIC devia…³
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.823 0.814 4.67 91.7 7.00e-29 4 -246. 504. 519. 1721.
## # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
## # variable names ¹adj.r.squared, ²statistic, ³deviance
broom::tidy(Fit1)
## # A tibble: 6 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -33.5 4.08 -8.22 3.22e-12
## 2 I(log(conc)) 8.48 0.678 12.5 2.04e-20
## 3 TypeQuebec:Treatmentnonchilled 19.5 1.44 13.6 2.59e-22
## 4 TypeMississippi:Treatmentnonchilled 10.1 1.44 7.04 6.25e-10
## 5 TypeQuebec:Treatmentchilled 15.9 1.44 11.1 9.91e-18
## 6 TypeMississippi:Treatmentchilled NA NA NA NA
broom.mixed::glance(Fit2)
## # A tibble: 1 × 7
## nobs sigma logLik AIC BIC REMLcrit df.residual
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 84 4.50 -241. 496. 513. 482. 77
broom.mixed::tidy(Fit2)
## # A tibble: 7 × 6
## effect group term estimate std.e…¹ stati…²
## <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) -33.5 4.02 -8.34
## 2 fixed <NA> I(log(conc)) 8.48 0.654 13.0
## 3 fixed <NA> TypeQuebec:Treatmentnonchilled 19.5 1.83 10.6
## 4 fixed <NA> TypeMississippi:Treatmentnonchilled 10.1 1.83 5.53
## 5 fixed <NA> TypeQuebec:Treatmentchilled 15.9 1.83 8.69
## 6 ran_pars Plant sd__(Intercept) 1.47 NA NA
## 7 ran_pars Residual sd__Observation 4.50 NA NA
## # … with abbreviated variable names ¹std.error, ²statistic
data("ChickWeight")
ggplot(ChickWeight, aes(x = Time, y = weight)) +
geom_point(aes(color = Diet)) +
geom_path(aes(color = Diet, group = Chick))
Fit1_Poisson <- glm(weight ~ Diet:Time, data =ChickWeight, family = poisson)
Fit2_Poisson <- glmer(weight ~ Diet:Time + (1|Chick), data =ChickWeight, family = poisson)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
broom::tidy(Fit1_Poisson)
## # A tibble: 5 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.86 0.00931 415. 0
## 2 Diet1:Time 0.0656 0.000722 90.8 0
## 3 Diet2:Time 0.0754 0.000768 98.1 0
## 4 Diet3:Time 0.0862 0.000729 118. 0
## 5 Diet4:Time 0.0823 0.000756 109. 0
broom.mixed::tidy(Fit2_Poisson)
## # A tibble: 6 × 7
## effect group term estimate std.error statistic p.value
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 3.84 0.0315 122. 0
## 2 fixed <NA> Diet1:Time 0.0666 0.00105 63.7 0
## 3 fixed <NA> Diet2:Time 0.0749 0.00129 57.9 0
## 4 fixed <NA> Diet3:Time 0.0868 0.00123 70.5 0
## 5 fixed <NA> Diet4:Time 0.0763 0.00126 60.4 0
## 6 ran_pars Chick sd__(Intercept) 0.213 NA NA NA
sjPlot::plot_model(Fit1_Poisson,type = "eff")
## Package `effects` is not available, but needed for `ggeffect()`. Either
## install package `effects`, or use `ggpredict()`. Calling `ggpredict()`
## now.
## $Diet
##
## $Time
para que sea un modelo aleatorio se debe tomar mas d euna medida por individuo