library('readxl')
df1<- read_excel(file.choose())
df <- df1

library('dplyr')  
df <- df %>%  
  mutate(  
    L6 = case_when(  
      L6 == "G" ~ 1,  
      L6 == "G/A" ~ 2,  
      L6 == "A" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),  
    L4 = case_when(  
      L4 == "T" ~ 1,  
      L4 == "C/T" ~ 2,  
      L4 == "C" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L3 = case_when(  
      L3 == "C" ~ 1,  
      L3 == "C/T" ~ 2,  
      L3 == "T" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L2 = case_when(  
      L2 == "C" ~ 1,  
      L2 == "C/G" ~ 2,  
      L2 == "G" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L1 = case_when(  
      L1 == "T" ~ 1,  
      L1 == "C/T" ~ 2,  
      L1 == "C" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L10 = case_when(  
      L10 == "C" ~ 1,  
      L10 == "C/T" ~ 2,  
      L10 == "T" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L9 = case_when(  
      L9 == "C" ~ 1,  
      L9 == "T/C" ~ 2,  
      L9 == "T" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L5 = case_when(  
      L5 == "A" ~ 1,  
      L5 == "G/A" ~ 2,  
      L5 == "G" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L8 = case_when(  
      L8 == "C" ~ 1,  
      L8 == "C/G" ~ 2,  
      L8 == "G" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L15 = case_when(  
      L15 == "C" ~ 1,  
      L15 == "C/T" ~ 2,  
      L15 == "T" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L13 = case_when(  
      L13 == "G" ~ 1,  
      L13 == "G/C" ~ 2,  
      L13== "C" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L12 = case_when(  
      L12 == "G" ~ 1,  
      L12 == "T/G" ~ 2,  
      L12== "T" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L17 = case_when(  
      L17 == "C" ~ 1,  
      L17 == "C/T" ~ 2,  
      L17== "T" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L14 = case_when(  
      L14 == "C" ~ 1,  
      L14 == "A/C" ~ 2,  
      L14 == "A" ~ 3,  
      TRUE ~ as.numeric(NA)  
    ),
    L16 = case_when(  
      L16 == "G" ~ 1,  
      L16 == "G/A" ~ 2,  
      L16== "A" ~ 3,  
      TRUE ~ as.numeric(NA)  
    )
  )

library(writexl)


df  <- df %>%  
  mutate_at(vars(-sample), as.factor)  
str(df)

df$Ph_1<- factor(df$Ph_1, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_2<- factor(df$Ph_2, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_3<- factor(df$Ph_3, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_4<- factor(df$Ph_4, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_5<- factor(df$Ph_5, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_6<- factor(df$Ph_6, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_7<- factor(df$Ph_7, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_8<- factor(df$Ph_8, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_9<- factor(df$Ph_9, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_10<- factor(df$Ph_10, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_11<- factor(df$Ph_11, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_12<- factor(df$Ph_12, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))
df$Ph_13<- factor(df$Ph_13, levels = c(2,1,0), labels = c("level_2", "level_1", "level_0"))


library('nnet')
model_multinom <- multinom(Similarity~L1 + L2 + L3 + L4+ L5 + L6 + L8 + L9 + L10 +   
                             L12 + L13 + L14 + L15 +   
                             L16 + L17,
                           data = df, maxit = 500)


aic_value <- model_multinom$AIC
print(paste("AIC值:", aic_value))

zvalues<-summary(model_multinom)$coefficients/
  summary(model_multinom)$standard.errors
p_mul <- pnorm(abs (zvalues), lower.tail = FALSE) * 2
p_mul_df <- as.data.frame(p_mul) 

coef_mul<- summary(model_multinom)$coefficients
#coef_mul_df <- as.data.frame(coef_mul) 

std.error_mul<- summary(model_multinom)$standard.errors
#std.error_mul_df <- as.data.frame(std.error_mul) 

OR_mul <- exp(coef_mul)
#OR_mul_df <- as.data.frame(OR_mul) 


ci <- confint(model_multinom )  
ci_mul_df <- as.data.frame(ci)