Module 15d

# Load required libraries
library(lme4)

Loading required package: Matrix

library(ggplot2)

# Check files and data structure
print(getwd())

[1] "/Users/seanlawler/Downloads/Intro to R/Module 15"

print(list.files())

 [1] "~$dule 15.docx"            "15a_files"                
 [3] "15a.html"                  "15a.qmd"                  
 [5] "15b_files"                 "15b.html"                 
 [7] "15b.qmd"                   "15c_files"                
 [9] "15c.html"                  "15c.qmd"                  
[11] "15d.qmd"                   "15d.rmarkdown"            
[13] "ca_places.zip"             "gapminderData5.csv"       
[15] "GEOG_5680_15a.html"        "GEOG_5680_15b.html"       
[17] "GEOG_5680_15c.html"        "GEOG_5680_15d.html"       
[19] "hsa.csv"                   "irished.csv"              
[21] "LightRail_UTA"             "LightRail_UTA.zip"        
[23] "LightRailStations_UTA"     "LightRailStations_UTA.zip"
[25] "Module 15.docx"            "Module 15.pdf"            
[27] "Module15_abc.qmd"          "Module15a_files"          
[29] "Module15a-b_files"         "Module15a-b.html"         
[31] "Module15a.html"            "orthodont.csv"            
[33] "rs.zip"                    "rsconnect"                
[35] "slc_tract"                 "slc_tract.zip"            
[37] "WNAclimate.csv"           

# Load and prepare data
ortho <- read.csv("orthodont.csv")
ortho$Subject <- factor(ortho$Subject)  # convert to factor

# Inspect structure
str(ortho)

'data.frame':   108 obs. of  4 variables:
 $ distance: num  26 25 29 31 21.5 22.5 23 26.5 23 22.5 ...
 $ age     : int  8 10 12 14 8 10 12 14 8 10 ...
 $ Subject : Factor w/ 27 levels "F01","F02","F03",..: 12 12 12 12 13 13 13 13 14 14 ...
 $ Sex     : chr  "Male" "Male" "Male" "Male" ...

# Model 1: Random Intercept Model
model1 <- lmer(distance ~ age + (1 | Subject), data = ortho, REML = FALSE)
summary(model1)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: distance ~ age + (1 | Subject)
   Data: ortho

      AIC       BIC    logLik -2*log(L)  df.resid 
    451.4     462.1    -221.7     443.4       104 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.6870 -0.5386 -0.0123  0.4910  3.7470 

Random effects:
 Groups   Name        Variance Std.Dev.
 Subject  (Intercept) 4.294    2.072   
 Residual             2.024    1.423   
Number of obs: 108, groups:  Subject, 27

Fixed effects:
            Estimate Std. Error t value
(Intercept) 16.76111    0.79456   21.09
age          0.66019    0.06122   10.78

Correlation of Fixed Effects:
    (Intr)
age -0.848

# Model 2: Random Slope Model
model2 <- lmer(distance ~ age + (age | Subject), data = ortho, REML = FALSE)
summary(model2)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: distance ~ age + (age | Subject)
   Data: ortho

      AIC       BIC    logLik -2*log(L)  df.resid 
    451.2     467.3    -219.6     439.2       102 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.3060 -0.4874  0.0076  0.4822  3.9228 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Subject  (Intercept) 4.81397  2.1941        
          age         0.04619  0.2149   -0.58
 Residual             1.71623  1.3100        
Number of obs: 108, groups:  Subject, 27

Fixed effects:
            Estimate Std. Error t value
(Intercept) 16.76111    0.76076  22.032
age          0.66019    0.06992   9.442

Correlation of Fixed Effects:
    (Intr)
age -0.848

# Model Comparison
anova(model1, model2)

Data: ortho
Models:
model1: distance ~ age + (1 | Subject)
model2: distance ~ age + (age | Subject)
       npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
model1    4 451.39 462.12 -221.69    443.39                     
model2    6 451.21 467.30 -219.61    439.21 4.1779  2     0.1238

# Individual growth curves
ggplot(ortho, aes(x = age, y = distance, group = Subject)) +
  geom_line(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE, color = "blue") +
  labs(title = "Growth of Dental Distance by Age",
       x = "Age (years)", y = "Distance (mm)")

`geom_smooth()` using formula = 'y ~ x'

---