# ============================================================

# এখানে CRD, RBD, Factorial Design, ANCOVA এবং Confounding ডিজাইনের বিশ্লেষণ রয়েছে

# ============================================================

## ==================== Factorial Design বিশ্লেষণ ====================
# Workspace ক্লিয়ার করা (সব ভ্যারিয়েবল ও ডাটা রিমুভ)
rm(list=ls())

# ডাটা লোড করা
Data1=read.csv("3^2 data.csv")
head(Data1)  # প্রথম কয়েক সারি দেখা

# dplyr লাইব্রেরি লোড করা
library(dplyr)

# Yield ছাড়া সব কলামকে ফ্যাক্টরে কনভার্ট করা
Data_RBD_Con=Data1 %>%
  mutate(across(!"Yield", as.factor))

# CRD মডেল ফিট করা (3-way factorial design)
model_CRD=aov(Yield~Nozzle_Type*Speed*Pressure, data = Data_RBD_Con)
summary(model_CRD)  # ANOVA টেবিল দেখা

# Nozzle_Type এর জন্য পোস্ট-হক টেস্ট (Tukey HSD)
TukeyHSD(model_CRD, "Nozzle_Type")
pairwise.t.test(Data_RBD_Con$Yield, Data_RBD_Con$Nozzle_Type, p.adjust.method =
                  "bonferroni")

# Speed এর জন্য পোস্ট-হক টেস্ট
TukeyHSD(model_CRD, "Speed")
pairwise.t.test(Data_RBD_Con$Yield, Data_RBD_Con$Speed, p.adjust.method =
                  "bonferroni")

# Pressure এর জন্য পোস্ট-হক টেস্ট
TukeyHSD(model_CRD, "Pressure")
pairwise.t.test(Data_RBD_Con$Yield, Data_RBD_Con$Pressure, p.adjust.method =
                  "bonferroni")

# গ্রাফিক্যাল ডিসপ্লের জন্য মার্জিন সেট করা
par(mar = c(5, 8, 4, 2)) # লেবেল স্পেসিং এর জন্য মার্জিন সেট

# Tukey HSD ফলাফল প্লট করা
plot(TukeyHSD(model_CRD, "Nozzle_Type"), las = 1, col = "pink")
plot(TukeyHSD(model_CRD, "Speed"), las = 1, col = "blue")
plot(TukeyHSD(model_CRD, "Pressure"), las = 1, col = "blue")

# Speed এর জন্য Mean Yield গ্রাফ
mean_data <- Data_RBD_Con %>%
  group_by(Speed) %>%
  summarise(mean_yield = mean(Yield))

attach(mean_data)
plot(mean_yield,Speed, data = mean_data, pch = 19, col = "blue",
     xlab = "Speed", ylab = "Mean Yield", main = "Mean Yield by Speed")

# আবার একই কোড (মূল কোডে দুবার আছে, তাই রিপিট করছি)
mean_data <- Data_RBD_Con %>%
  group_by(Speed) %>%
  summarise(mean_yield = mean(Yield))

attach(mean_data)
plot(mean_yield,Speed, data = mean_data, pch = 19, col = "blue",
     xlab = "Speed", ylab = "Mean Yield", main = "Mean Yield by Speed")


## ==================== ANCOVA (Analysis of Covariance) ====================
# Workspace ক্লিয়ার করা
rm(list=ls())

# কারেন্ট ওয়ার্কিং ডিরেক্টরি চেক করা
getwd()

# ANCOVA ডাটা লোড করা
data<-data.frame(read.csv("ancova.csv"))
data  # ডাটা দেখানো

# Potato ভ্যারিয়েবলকে ফ্যাক্টরে কনভার্ট করা
data$Potatos=as.factor(data$Potato)

# ANOVA মডেল ফিট করা (X এবং Yield এর জন্য আলাদা)
av1=aov(X~Potatos,data=data)   # কভ্যারিয়েটের জন্য
av2=aov(Yield~Potatos,data=data)  # রেসপন্স ভ্যারিয়েবলের জন্য

# প্রথম ANOVA সামারি
summary(av1)

# ANOVA টেবিল থেকে Sum of Squares বের করা
Txx=anova(av1)["Potatos","Sum Sq"]  # ট্রিটমেন্ট সমষ্টি বর্গ
Exx=anova(av1)["Residuals","Sum Sq"]  # এরর সমষ্টি বর্গ
Gxx=Txx+Exx  # টোটাল সমষ্টি বর্গ

Tyy=anova(av2)["Potatos","Sum Sq"]  # Yield এর ট্রিটমেন্ট সমষ্টি বর্গ
Eyy=anova(av2)["Residuals","Sum Sq"]  # Yield এর এরর সমষ্টি বর্গ

# দ্বিতীয় ANOVA সামারি
summary(av2)

Gyy=Tyy+Eyy  # Yield এর টোটাল সমষ্টি বর্গ

# রিগ্রেশন মডেল ফিট করা
reg1=lm(Yield~Potatos+X,data=data)  # ANCOVA মডেল
reg2=lm(Yield~X,data=data)  # সরল রিগ্রেশন মডেল

# রিগ্রেশন কোএফিশিয়েন্ট বের করা
b_hat=coef(reg1)["X"]  # প্রথম উপায়
b_hat=coef(reg1)[5]    # ভেক্টর ইন্ডেক্স দিয়ে (দ্বিতীয় উপায়)
b_2_hat=coef(reg2)[2]  # দ্বিতীয় মডেলের কোএফিশিয়েন্ট

# Cross-products বের করা
Exy=Exx*b_hat   # এরর ক্রস প্রোডাক্ট
Gxy=Gxx*b_2_hat  # টোটাল ক্রস প্রোডাক্ট
Txy=Gxy-Exy     # ট্রিটমেন্ট ক্রস প্রোডাক্ট

# ANCOVA মডেল ফিট করা
ancova_model <- aov(Yield ~ Potato + X, data = data)
summary(ancova_model)  # ANCOVA ফলাফল দেখা

# ---- ধাপ ৬: রিগ্রেশন ঢালের homogeneity চেক করা ----
# ইন্টার‍্যাকশন টার্ম সহ মডেল (ঢালের homogeneity টেস্টের জন্য)
ancova_interaction <- aov(Yield ~ Potato * X, data = data)
summary(ancova_interaction)

# মডেল কম্পেয়ার করা
anova(ancova_model, ancova_interaction)

# ---- ধাপ ৭: ডায়াগনস্টিক প্লট ----
par(mfrow = c(2, 2))  # ২x২ গ্রিডে প্লট সেট করা
plot(ancova_model)    # রেসিডুয়াল প্লট

# emmeans লাইব্রেরি লোড করা (পোস্ট-হক টেস্টের জন্য)
library(emmeans)
emmeans(ancova_model, "Potato")  # adjusted means দেখা


## ==================== Confounding ডিজাইন বিশ্লেষণ ====================
# কনফাউন্ডিং ডিজাইনের ডাটা লোড করা
data<-read.csv("onnesha_confounding.csv")
data  # ডাটা দেখা

# সব ভ্যারিয়েবলকে (Yield ছাড়া) ফ্যাক্টরে কনভার্ট করা
data_rbd= data %>%
  mutate(across(! "Yield",as.factor))

# RBD মডেল ফিট করা (কনফাউন্ডিং সহ)
model_rbd=aov(Yield~Block+N*P*K*D,data=data_rbd)
summary(model_rbd)  # ANOVA টেবিল দেখা

# প্রতিটি ফ্যাক্টরের জন্য পোস্ট-হক টেস্ট
TukeyHSD(model_rbd,"N")
pairwise.t.test(data_rbd$Yield,data_rbd$N, p.adjust="bonferroni")

TukeyHSD(model_rbd,"P")
pairwise.t.test(data_rbd$Yield,data_rbd$P, p.adjust="bonferroni")

TukeyHSD(model_rbd,"K")
pairwise.t.test(data_rbd$Yield,data_rbd$K, p.adjust="bonferroni")

TukeyHSD(model_rbd,"D")
pairwise.t.test(data_rbd$Yield,data_rbd$D, p.adjust="bonferroni")

# Tukey HSD ফলাফল প্লট করা
plot(TukeyHSD(model_rbd,"N"), las=1, col="blue")
plot(TukeyHSD(model_rbd,"P"), las=1, col="red")
plot(TukeyHSD(model_rbd,"K"), las=1, col="blue")
plot(TukeyHSD(model_rbd,"D"), las=1, col="blue")

# N (নাইট্রোজেন) ফ্যাক্টরের জন্য Mean Yield গ্রাফ
mean_data<-data_rbd %>%
  group_by(N) %>%
  summarise(mean_yield=mean(Yield))

attach(mean_data)
plot(mean_yield, N, data=mean_data, pch=1, X_lab="Nitrogen", y_lab="Mean Yield")

# P (পটাশ) ফ্যাক্টরের জন্য Mean Yield গ্রাফ
mean_data<-data_rbd %>%
  group_by(P) %>%
  summarise(mean_yield=mean(Yield))

attach(mean_data)
plot(mean_yield, P, data=mean_data, pch=1, X_lab="Potash", y_lab="Mean Yield")


## ==================== ল্যাব ১ম ও ২য় ক্লাস ====================
## ==================== CRD (Completely Randomized Design) ====================
# Workspace ক্লিয়ার করা
rm(list = ls())

# CRD ডাটা লোড করা
data <- read.csv("CRD_404.csv")

# ডাটার প্রথম কয়েক সারি দেখা
head(data)

# ট্রিটমেন্ট ভ্যারিয়েবলগুলোকে ফ্যাক্টরে কনভার্ট করা
data$Treatment_A <- as.factor(data$Treatment_A)
data$Treatment_B <- as.factor(data$Treatment_B)

# টু-ওয়ে ANOVA with ইন্টার‍্যাকশন ফিট করা
model <- aov(Yield ~ Treatment_A * Treatment_B, data = data)

# ANOVA টেবিল দেখানো
summary(model)

# Type II বা Type III Sum of Squares (Python এর statsmodels এর typ=2 এর মতো)
# car প্যাকেজ ইনস্টল করা লাগতে পারে (যদি না থাকে)
# install.packages("car")
library(car)

# Type II ANOVA (statsmodels এর typ=2 এর সমতুল্য)
Anova(model, type = 2)

# Type III ANOVA (যদি ডিজাইন আনব্যালান্সড হয়)
Anova(model, type = 3)

# ঐচ্ছিক: ইন্টার‍্যাকশন প্লট দেখা
interaction.plot(data$Treatment_A, data$Treatment_B, data$Yield,
                 main = "Interaction Plot",
                 xlab = "Treatment A", ylab = "Yield", col = 2:4, lwd = 2)

# ঐচ্ছিক: ডায়াগনস্টিক প্লট
par(mfrow = c(2, 2))
plot(model)


## ==================== RBD (Randomized Block Design) ====================
# প্রথম RBD উদাহরণ
data<-read.csv("RBD 1.csv")
data  # ডাটা দেখা

# ভ্যারিয়েবলগুলোকে ফ্যাক্টরে কনভার্ট করা
block=as.factor(data$Block)
age_group= as.factor(data$Age_group)
body_condition=as.factor(data$Body_condition)
protein_level= as.factor(data$Protein)

# RBD মডেল ফিট করা (ব্লক সহ)
Model_RBD=aov(Body_weight~ block+age_group*body_condition*protein_level,
              data=data)
summary(Model_RBD)  # ANOVA ফলাফল দেখা


## ==================== RBD ২য় উদাহরণ ====================
# Workspace ক্লিয়ার করা
rm(list = ls())

# RBD ডাটা লোড করা
data <- read.csv("RBD_404 - Copy.csv")
getwd()  # কারেন্ট ওয়ার্কিং ডিরেক্টরি চেক

# ক্যাটেগরিক্যাল ভ্যারিয়েবলগুলোকে ফ্যাক্টরে কনভার্ট করা
data$Block <- as.factor(data$Block)
data$Treatment_A <- as.factor(data$Treatment_A)
data$Treatment_B <- as.factor(data$Treatment_B)

# ব্লকিং ও ইন্টার‍্যাকশন সহ টু-ওয়ে ANOVA মডেল ফিট করা
model <- aov(Yield ~ Block + Treatment_A * Treatment_B, data = data)

# ANOVA টেবিল দেখানো (ডিফল্ট Type I SS)
summary(model)

# Type II SS এর জন্য (Python এর statsmodels এর typ=2 এর মতো)
# install.packages("car") # একবার রান করলেই হবে যদি না ইন্সটল থাকে
library(car)
Anova(model, type = 2)

# ঐচ্ছিক: ডায়াগনস্টিক প্লট
par(mfrow = c(2, 2))
plot(model)