Running examples in 'DImodels-Ex.R' failed
The error most likely occurred in:

> ### Name: sim2
> ### Title: The Simulated "sim2" Dataset
> ### Aliases: sim2
> ### Keywords: datasets
> 
> ### ** Examples
> 
> 
> ####################################
> ## Code to simulate the sim2 dataset
>   
> 
> 
> 
> ###################################################################################################
> ###################################################################################################
> ## sim2
> 
> ###########################
> ## Analyse the sim2 dataset
>   
> ## Load the sim2 data
>   data(sim2)
> ## View the first few entries
>   head(sim2)
  community block  p1  p2  p3  p4 response
1         1     1 0.7 0.1 0.1 0.1   20.230
2         1     2 0.7 0.1 0.1 0.1   20.070
3         1     3 0.7 0.1 0.1 0.1   20.858
4         1     4 0.7 0.1 0.1 0.1   17.044
5         2     1 0.1 0.7 0.1 0.1   17.159
6         2     2 0.1 0.7 0.1 0.1   19.899
> ## Explore the variables in sim2
>   str(sim2)
'data.frame':	60 obs. of  7 variables:
 $ community: int  1 1 1 1 2 2 2 2 3 3 ...
 $ block    : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
 $ p1       : num  0.7 0.7 0.7 0.7 0.1 0.1 0.1 0.1 0.1 0.1 ...
 $ p2       : num  0.1 0.1 0.1 0.1 0.7 0.7 0.7 0.7 0.1 0.1 ...
 $ p3       : num  0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.7 0.7 ...
 $ p4       : num  0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
 $ response : num  20.2 20.1 20.9 17 17.2 ...
>   
> ## Check that the proportions sum to 1 (required for DI models)
> ## p1 to p4 are in the 3rd to 6th columns in sim2
>   sim2sums <- rowSums(sim2[3:6])
>   summary(sim2sums)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1       1       1       1       1       1 
>   
> ## Check characteristics of sim2
>   hist(sim2$response)
>   summary(sim2$response)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   6.36   12.03   18.37   16.67   20.27   22.45 
>   plot(sim2$p1, sim2$response)
>   plot(sim2$p2, sim2$response)
>   plot(sim2$p3, sim2$response)
>   plot(sim2$p4, sim2$response)
> 
> ## Find the best DI model using autoDI and F-test selection
>   auto1 <- autoDI(y = "response", prop = c("p1", "p2", "p3", "p4"), block = "block", data = sim2, 
+                   selection = "Ftest")

--------------------------------------------------------------------------------
Step 1: Investigating the diversity effect

All models include block

Selection using F tests

           Description                                                 
DI Model 1 Structural 'STR' DImodel                                    
DI Model 2 Species identity 'ID' DImodel                               
DI Model 3 Average interactions 'AV' DImodel                           
DI Model 4 Additive species contributions to interactions 'ADD' DImodel
DI Model 5 Separate pairwise interactions 'FULL' DImodel               


           DI_model treat estimate_theta Resid. Df Resid. SSq Resid. MSq Df
DI Model 1      STR  none          FALSE        56  1296.9569    23.1599   
DI Model 2       ID  none          FALSE        53  1267.5762    23.9165  3
DI Model 3       AV  none          FALSE        52    73.9632     1.4224  1
DI Model 4      ADD  none          FALSE        49    70.0529     1.4297  3
DI Model 5     FULL  none          FALSE        47    68.8948     1.4658  2
                SSq        F  Pr(>F)
DI Model 1                          
DI Model 2  29.3807   6.6812   8e-04
DI Model 3 1193.613 814.2827 <0.0001
DI Model 4   3.9103   0.8892  0.4537
DI Model 5   1.1581    0.395  0.6759

Functional groups (argument 'FG') were not specified, and therefore not investigated.

Selected model: Average interactions 'AV' DImodel
Formula: 
response  =  p1 + p2 + p3 + p4 + AV + block1 + block2 + block3 

--------------------------------------------------------------------------------
Step 2: No investigation of treatment effect included, since no treatment was specified
        (argument 'treat' omitted)

--------------------------------------------------------------------------------
Step 3: Investigating whether theta is equal to 1 or not

Theta estimate: 0.4533

Selection using F tests

           Description                                        
DI Model 1 Average interactions 'AV' DImodel                  
DI Model 2 Average interactions 'AV' DImodel, estimating theta


           DI_model treat estimate_theta Resid. Df Resid. SSq Resid. MSq Df
DI Model 1       AV  none          FALSE        52    73.9632     1.4224   
DI Model 2       AV  none           TRUE        51    55.0925     1.0802  1
               SSq       F Pr(>F)
DI Model 1                       
DI Model 2 18.8707 17.4689  1e-04

Selected model: Average interactions 'AV' DImodel, estimating theta


--------------------------------------------------------------------------------
Step 4: Comparing the final selected model with the reference (community) model
'community' is a factor with 15 levels, one for each unique set of proportions.



               model Resid. Df Resid. SSq Resid. MSq Df   SSq      F Pr(>F)
DI Model 1  Selected        51    55.0925     1.0802                       
DI Model 2 Reference        42    48.8616     1.1634  9 6.231 0.5951 0.7936

--------------------------------------------------------------------------------
autoDI is limited in terms of model selection. Exercise caution when choosing your final model.
--------------------------------------------------------------------------------
>   summary(auto1)

Call:
glm(formula = formula(obj), family = family, data = data_theta_E_AV)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.70854  -0.74378   0.05998   0.67862   2.48735  

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
p1       9.3837     0.5015  18.713  < 2e-16 ***
p2       8.3691     0.5015  16.690  < 2e-16 ***
p3       6.9869     0.5015  13.933  < 2e-16 ***
p4       6.6656     0.5015  13.293  < 2e-16 ***
AV       6.9954     0.2088  33.502  < 2e-16 ***
block1   1.4945     0.3795   3.938  0.00025 ***
block2   2.2542     0.3795   5.940 2.55e-07 ***
block3   2.1809     0.3795   5.746 5.11e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.080246)

    Null deviance: 18014.041  on 60  degrees of freedom
Residual deviance:    55.093  on 51  degrees of freedom
AIC: 185.15

Number of Fisher Scoring iterations: 2

>   
>   ## Fit the average pairwise model, including theta, using DI and the AV tag
>   m1 <- DI(y = "response", prop = c("p1","p2","p3","p4"), block = "block", DImodel = "AV", 
+            estimate_theta = TRUE, data = sim2)
Fitted model: Average interactions 'AV' DImodel
Theta estimate: 0.4533
>   summary(m1)

Call:
glm(formula = formula(obj), family = family, data = data_theta_E_AV)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.70854  -0.74378   0.05998   0.67862   2.48735  

Coefficients:
       Estimate Std. Error t value Pr(>|t|)    
p1       9.3837     0.5015  18.713  < 2e-16 ***
p2       8.3691     0.5015  16.690  < 2e-16 ***
p3       6.9869     0.5015  13.933  < 2e-16 ***
p4       6.6656     0.5015  13.293  < 2e-16 ***
AV       6.9954     0.2088  33.502  < 2e-16 ***
block1   1.4945     0.3795   3.938  0.00025 ***
block2   2.2542     0.3795   5.940 2.55e-07 ***
block3   2.1809     0.3795   5.746 5.11e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 1.080246)

    Null deviance: 18014.041  on 60  degrees of freedom
Residual deviance:    55.093  on 51  degrees of freedom
AIC: 185.15

Number of Fisher Scoring iterations: 2

>   CI_95 <- theta_CI(m1, conf = .95)
>   CI_95
    lower     upper 
0.2861322 0.6674397 
>   plot(m1)
>   library(hnp)
Loading required package: MASS
> 
> 
> ## Graph the profile likelihood
>   library(ggplot2)
Error in library(ggplot2) : there is no package called 'ggplot2'
Execution halted