### Practicals 11 
library(ggplot2)
library(readxl)

### task 1 
grass<-read_excel("11_data.xlsx", sheet="1")
summary(grass)

#These are observatory data; therefore, it is reasonable to use the forward-selection
# approach
lm.0<-lm(seedlings~+1, data=grass)# Null (intercept-only) model

add1(lm.0, .~.+treatment*productivity*temperature, test="F")# testing addition of individual candidate
# predictors. Interactions are also included here but will only be tested when both main effects of
#corresponding interaction pass in.
# 
# Single term additions
# 
# Model:
# seedlings ~ +1
#              Df Sum of Sq    RSS     AIC F value    Pr(>F)    
# <none>                    819.37 101.220                      
# treatment     2    400.47 418.90  85.093 12.9059 0.0001166 ***
# productivity  1     81.11 738.25 100.093  3.0765 0.0903746 .  
# temperature   1     39.55 779.81 101.736  1.4203 0.2433658    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lm.1<-update(lm.0, .~.+treatment)# The best predictor (lowest AIC) is treatment and thus is added to the model
# Note this part of the code: ".~.+treatment" it meens leave the model as is just add the treatment predictor
add1(lm.1, .~.+treatment*productivity*temperature, test="F")
# Single term additions
# 
# Model:
# seedlings ~ treatment
#              Df Sum of Sq    RSS    AIC F value   Pr(>F)   
# <none>                    418.90 85.093                    
# productivity  1    99.658 319.24 78.943  8.1165 0.008465 **
# temperature   1    45.090 373.81 83.677  3.1362 0.088297 . 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lm.2<-update(lm.1, .~.+productivity)
add1(lm.2, .~.+treatment*productivity*temperature, test="F")
# Single term additions
# 
# Model:
# seedlings ~ treatment + productivity
#                        Df Sum of Sq    RSS    AIC F value Pr(>F)  
# <none>                              319.24 78.943                 
# temperature             1    25.781 293.46 78.416  2.1963 0.1508  
# treatment:productivity  2    69.397 249.84 75.589  3.3331 0.0528 .
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(lm.2)# Checking model assumptions; Seems okay

anova(lm.2)
# Analysis of Variance Table
# 
# Response: seedlings
#              Df Sum Sq Mean Sq F value    Pr(>F)    
# treatment     2 400.47 200.233 16.3076 2.573e-05 ***
# productivity  1  99.66  99.658  8.1165  0.008465 ** 
# Residuals    26 319.24  12.279                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


summary(lm.2)
# Analysis of Variance Table
# 
# Response: seedlings
#              Df Sum Sq Mean Sq F value    Pr(>F)    
# treatment     2 400.47 200.233 16.3076 2.573e-05 ***
# productivity  1  99.66  99.658  8.1165  0.008465 ** 
# Residuals    26 319.24  12.279                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# > summary(lm.2)
# 
# Call:
# lm(formula = seedlings ~ treatment + productivity, data = grass)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -7.0899 -2.6108  0.0938  1.7933  6.8140 
# 
# Coefficients:
#                 Estimate Std. Error t value Pr(>|t|)    
# (Intercept)    25.105812   2.538236   9.891 2.66e-10 ***
# treatmentmown1 -7.882879   1.568384  -5.026 3.14e-05 ***
# treatmentmown2 -7.990293   1.568492  -5.094 2.62e-05 ***
# productivity   -0.012357   0.004337  -2.849  0.00847 ** 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 3.504 on 26 degrees of freedom
# Multiple R-squared:  0.6104,	Adjusted R-squared:  0.5654 
# F-statistic: 13.58 on 3 and 26 DF,  p-value: 1.594e-05



#Task 2
transpi<-read_excel("11_data.xlsx", sheet="2")
summary(transpi)


#Here, again we work with observatory data; therefore forward selection is a reasonable approach. 

lm.0<-lm(transpiration~+1, data=transpi)
add1(lm.0, ~+wind*temperature*humidity*sunshine, test="F")
# Single term additions
# 
# Model:
# transpiration ~ +1
#             Df Sum of Sq     RSS    AIC F value   Pr(>F)   
# <none>                   158.000 35.930                    
# wind         1    75.369  82.631 28.855 10.9454 0.006243 **
# temperature  1    93.763  64.237 25.329 17.5157 0.001265 **
# humidity     1     0.038 157.962 37.926  0.0029 0.958179   
# sunshine     1     4.667 153.333 37.510  0.3652 0.556876   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
lm.1<-update(lm.0, .~.+temperature)
add1(lm.1, ~+wind*temperature*humidity*sunshine, test="F")
# Single term additions
# 
# Model:
# transpiration ~ temperature
#          Df Sum of Sq    RSS    AIC F value   Pr(>F)   
# <none>                64.237 25.329                    
# wind      1    37.523 26.714 15.046 15.4509 0.002349 **
# humidity  1     0.002 64.235 27.329  0.0004 0.984248   
# sunshine  1     0.229 64.009 27.279  0.0393 0.846520   
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
lm.2<-update(lm.1, .~.+wind)
add1(lm.2, ~+wind*temperature*humidity*sunshine, test="F")
# Single term additions
# 
# Model:
# transpiration ~ temperature + wind
#                  Df Sum of Sq    RSS    AIC F value Pr(>F)
# <none>                        26.714 15.046               
# humidity          1   2.60941 24.105 15.607  1.0825 0.3226
# sunshine          1   0.59063 26.123 16.733  0.2261 0.6447
# wind:temperature  1   0.81632 25.898 16.611  0.3152 0.5869

plot(lm.2)# Checking the residuals

summary(lm.2)
# Call:
# lm(formula = transpiration ~ temperature + wind, data = transpi)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -2.8698 -0.7002 -0.2190  0.8018  2.2732 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  5.45464    1.14007   4.784 0.000567 ***
# temperature  0.38619    0.08048   4.798 0.000555 ***
# wind         0.67796    0.17247   3.931 0.002349 ** 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 1.558 on 11 degrees of freedom
# Multiple R-squared:  0.8309,	Adjusted R-squared:  0.8002 
# F-statistic: 27.03 on 2 and 11 DF,  p-value: 5.681e-05

## Conclusion: There are significant positive effects of temperature and wind on transpiration.

# Task 3

plants<-read_excel("11_data.xlsx", sheet="3")
summary(plants)

# These are experimental data, so it make sense to fit a saturated model with interactions first
# and then remove possible non/significant predictors
lm.f<-lm(biomass~watering*fertil, data=plants)
anova(lm.f)
# Analysis of Variance Table
# 
# Response: biomass
#                 Df  Sum Sq Mean Sq F value    Pr(>F)    
# watering         1  6144.0  6144.0  23.345 0.0001013 ***
# fertil           1 14113.5 14113.5  53.626 4.443e-07 ***
# watering:fertil  1  4108.2  4108.2  15.610 0.0007891 ***
# Residuals       20  5263.7   263.2                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Both experimental treatments and their interaction are significant. There is no need to remove any
# predictor from the model

# This is also obvious here from the output of the function drop1  
drop1(lm.f, test="F")
# Single term deletions
# 
# Model:
# biomass ~ watering * fertil
#                 Df Sum of Sq    RSS    AIC F value    Pr(>F)    
# <none>                       5263.7 137.37                      
# watering:fertil  1    4108.2 9371.8 149.22   15.61 0.0007891 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Checking the assumptions
plot(lm.f)# There seems to be positive correlation between absolute values of residuals and
# predicted values, therefore, it may be better to refit the model with log-transformed data


lm.f2<-lm(log(biomass)~watering*fertil, data=plants)
anova(lm.f2)
# Analysis of Variance Table
# 
# Response: log(biomass)
#                 Df  Sum Sq Mean Sq F value    Pr(>F)    
# watering         1 0.41320 0.41320  24.730 7.314e-05 ***
# fertil           1 1.12613 1.12613  67.399 7.797e-08 ***
# watering:fertil  1 0.20946 0.20946  12.537  0.002052 ** 
# Residuals       20 0.33417 0.01671                      
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(lm.f2) # The residuals look much better now

summary(lm.f2)
# lm(formula = log(biomass) ~ watering * fertil, data = plants)
# 
# Residuals:
#       Min        1Q    Median        3Q       Max 
# -0.215320 -0.065523  0.009066  0.081658  0.213493 
# 
# Coefficients:
#                                  Estimate Std. Error t value Pr(>|t|)    
# (Intercept)                       4.36385    0.05277  82.695  < 2e-16 ***
# wateringWatered                   0.07558    0.07463   1.013  0.32328    
# fertilFertilized                  0.24639    0.07463   3.301  0.00356 ** 
# wateringWatered:fertilFertilized  0.37369    0.10554   3.541  0.00205 ** 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.1293 on 20 degrees of freedom
# Multiple R-squared:  0.8396,	Adjusted R-squared:  0.8155 
# F-statistic: 34.89 on 3 and 20 DF,  p-value: 3.862e-08

# Conclusion: Plant growth is significantly positively affect by fertilizer and its interaction with increased irrigation.
# Biologically speaking, this means that fertilizer supports growth in all conditions but 
# increased watering only if fertilizer is applied. 

#Plotting the interaction plot
ggplot(plants, aes(x=watering, y=biomass, col=fertil))+
  stat_summary(fun.data = mean_cl_normal)+scale_y_log10()+
  labs(col="Fertilizer application")+ theme_bw()


# Task 4

# Try yourself! :-) 


