Rent

rm(list = ls())
library(brms)
library(dplyr)
library(ggplot2)
library(pammtools)
library(grid)
library(gridExtra)
library(tidybayes)
library(bayesplot)

data("rent99", package = "gamlss.data")
head(rent99)
skimr::skim(rent99)

ggplot(rent99, aes(x=yearc, y=area, fill=rent))+
  geom_tile()+
  theme_bw()

ggplot(rent99, aes(x=yearc, 
                   y=area, 
                   size=rent,
                   color=as.factor(location)))+
  geom_point(alpha=0.5)+
  labs(color="Location")+
  theme_bw()
ggplot(rent99, aes(x=yearc, 
                   y=area, 
                   size=rent,
                   color=as.factor(location)))+
  geom_point(alpha=0.5)+
  labs(color="Location")+
  facet_wrap(~cheating)+
  theme_bw()
ggplot(rent99, aes(x=area, 
                   y=yearc, 
                   z=rentsqm))+
  geom_contour()+
  theme_bw()
ggplot(rent99, aes(x=area, 
                   y=yearc, 
                   z=rentsqm))+
  geom_density_2d()+
  facet_wrap(~cheating)+
  theme_bw()
ggplot(rent99, aes(x=area, 
                   y=yearc, 
                   z=rentsqm))+
  geom_density_2d_filled()+
  facet_wrap(~cheating)+
  theme_bw()


g <- mgcv::gam(rentsqm ~ te(area, yearc), 
               data=rent99)
plot(g)
gg_tensor(g)
gg_tensor(g, ci=TRUE)

get_prior(rentsqm ~ t2(area, yearc) + 
            (1|district), data = rent99, family=gaussian())
summary(rent99$rentsqm)
fit_rent1 <- brm(rentsqm ~ t2(area, yearc) + 
                   (1|district), 
                 data = rent99,
                 prior=c(
                 prior(normal(5, 2),  class = Intercept),
                 prior(gamma(10, 2),  class = sd,  coef= Intercept, group=district ),
                 prior(gamma(4, 2),  class = sds, coef= t2(area, yearc)),
                 prior(normal(3, 3),  class = b, coef=t2areayearc_1),
                 prior(normal(4, 4),  class = b, coef=t2areayearc_2),
                 prior(normal(3, 8),  class = b, coef=t2areayearc_3),
                 prior(gamma(5, 5),    class = sigma)),
                 chains = 4, cores = 6)
fit_rent1$prior
summary(fit_rent1)
plot(fit_rent1)
pairs(fit_rent1)
pp_check(fit_rent1)
ce<-conditional_effects(fit_rent1, surface = TRUE)
pl_ce<-plot(ce, rug=TRUE, theme=theme_bw())
grid.arrange(pl_ce$area,
             pl_ce$yearc, 
             pl_ce$`area:yearc`, 
             ncol=3)

bform <- bf(rentsqm ~ t2(area, yearc) + (1|ID1|district) +
              ar(p=3),
            sigma ~ t2(area, yearc) + (1|ID1|district))
get_prior(bform, data=rent99, family=gaussian())

g2 <- gamm4::gamm4(rentsqm ~ t2(area, yearc), 
                   random=~(1|district), 
               data=rent99)
g2
summary(g2$mer)
lkj_prior<-vcov(g2$mer)


fit_rent2 <- brm(bform, 
                 data = rent99,
                 prior=c(
                   prior(beta(3, 2),    class = ar),
                   prior(normal(3, 3),  class = b, coef=t2areayearc_1),
                   prior(normal(4, 4),  class = b, coef=t2areayearc_2),
                   prior(normal(3, 8),  class = b, coef=t2areayearc_3),
                   prior(lkj(1),  class = cor, group=district),
                   prior(normal(5, 2),  class = Intercept),
                   prior(gamma(10, 2),  class = sd,  coef= Intercept, group=district ),
                   prior(gamma(4, 2),   class = sds, coef= t2(area, yearc)),
                   prior(normal(3, 3),  class = b, coef=t2areayearc_1, dpar=sigma),
                   prior(normal(4, 4),  class = b, coef=t2areayearc_2, dpar=sigma),
                   prior(normal(3, 8),  class = b, coef=t2areayearc_3, dpar=sigma),
                   prior(gamma(10, 2),  class = sd,  coef= Intercept, group=district, dpar=sigma ),
                   prior(gamma(4, 2),   class = sds, coef= t2(area, yearc), dpar=sigma)),
                 chains = 4, cores = 6)
fit_rent2$prior
summary(fit_rent1)
summary(fit_rent2)
plot(fit_rent2)
pp_check(fit_rent2, ndraws=100)
loo_model1<-loo(fit_rent1)
loo_model2<-loo(fit_rent2)
loo_model1
loo_model2
loo_compare(loo_model1, loo_model2)
cs<-conditional_smooths(fit_rent2)
pl_cs<-plot(cs, theme=theme_bw())
grid.arrange(pl_cs$`mu: t2(area,yearc)`, 
             pl_cs$`sigma: t2(area,yearc)`, 
             ncol=2)
ce<-conditional_effects(fit_rent2, surface = TRUE)
pl_ce<-plot(ce, rug=TRUE, theme=theme_bw())