I will be analyzing a data frame containing data on The Consumption of Ice Cream. The variables include:

1.consumption of ice cream per head in pints 2.income of the average family per week 3.the price of ice cream per pint 4.and the average temperature in the U.S.

The observation numbers were 30 and it was across 4 weeks .

library(Ecdat)
library(Ecfun)
data("Icecream")
head(Icecream)
library(ggplot2)
ggplot(Icecream, aes(x=cons)) + geom_bar(fill="blue", width=.004)

So it seems as though the highest consumption of ice cream is about 2 pints.
ggplot(Icecream, aes(x=cons, y=temp)) + geom_bar(stat="identity",fill="blue", width=.004)

And then we can see here that as the temperature rises, so does the consumption of ice cream.

m1<-lm (cons~ price, data = Icecream)
summary(m1)

Call:
lm(formula = cons ~ price, data = Icecream)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.11152 -0.03707 -0.00991  0.03666  0.15724 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.9230     0.3964   2.329   0.0273 *
price        -2.0472     1.4393  -1.422   0.1660  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06466 on 28 degrees of freedom
Multiple R-squared:  0.06739,   Adjusted R-squared:  0.03408 
F-statistic: 2.023 on 1 and 28 DF,  p-value: 0.166
m2<-lm (cons ~ price + income, data = Icecream)
summary(m2)

Call:
lm(formula = cons ~ price + income, data = Icecream)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.10951 -0.03865 -0.01045  0.03665  0.15635 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept)  0.9002400  0.4550344   1.978   0.0582 .
price       -2.0300382  1.4738940  -1.377   0.1797  
income       0.0002135  0.0019687   0.108   0.9144  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06583 on 27 degrees of freedom
Multiple R-squared:  0.0678,    Adjusted R-squared:  -0.001257 
F-statistic: 0.9818 on 2 and 27 DF,  p-value: 0.3876
m3 <-lm (cons ~ price + temp, data = Icecream)
summary(m3)

Call:
lm(formula = cons ~ price + temp, data = Icecream)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.08226 -0.02051  0.00184  0.02272  0.10076 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.59655    0.25831   2.309   0.0288 *  
price       -1.40176    0.92509  -1.515   0.1413    
temp         0.00303    0.00047   6.448 6.56e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04132 on 27 degrees of freedom
Multiple R-squared:  0.6328,    Adjusted R-squared:  0.6056 
F-statistic: 23.27 on 2 and 27 DF,  p-value: 1.336e-06

Temperature has a strong significance on consumption in this linear model.

library(texreg)
screenreg(list(m1,m2,m3))

========================================
             Model 1  Model 2  Model 3  
----------------------------------------
(Intercept)   0.92 *   0.90     0.60 *  
             (0.40)   (0.46)   (0.26)   
price        -2.05    -2.03    -1.40    
             (1.44)   (1.47)   (0.93)   
income                 0.00             
                      (0.00)            
temp                            0.00 ***
                               (0.00)   
----------------------------------------
R^2           0.07     0.07     0.63    
Adj. R^2      0.03    -0.00     0.61    
Num. obs.    30       30       30       
RMSE          0.06     0.07     0.04    
========================================
*** p < 0.001, ** p < 0.01, * p < 0.05
m4<-lm(cons ~ price*temp, data=Icecream)
summary(m4)

Call:
lm(formula = cons ~ price * temp, data = Icecream)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.093826 -0.021935  0.006617  0.020962  0.082094 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.17029    0.89149  -0.191    0.850
price        1.38220    3.23280   0.428    0.672
temp         0.01880    0.01755   1.071    0.294
price:temp  -0.05731    0.06375  -0.899    0.377

Residual standard error: 0.04146 on 26 degrees of freedom
Multiple R-squared:  0.6439,    Adjusted R-squared:  0.6028 
F-statistic: 15.67 on 3 and 26 DF,  p-value: 5.07e-06
m0 <- glm(cons~ price, data = Icecream)
summary(m0)

Call:
glm(formula = cons ~ price, data = Icecream)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.11152  -0.03707  -0.00991   0.03666   0.15724  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   0.9230     0.3964   2.329   0.0273 *
price        -2.0472     1.4393  -1.422   0.1660  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.004180874)

    Null deviance: 0.12552  on 29  degrees of freedom
Residual deviance: 0.11706  on 28  degrees of freedom
AIC: -75.251

Number of Fisher Scoring iterations: 2
m01 <- glm(cons~ price + temp + income, data = Icecream)
summary(m01)

Call:
glm(formula = cons ~ price + temp + income, data = Icecream)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.065302  -0.011873   0.002737   0.015953   0.078986  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.1973151  0.2702162   0.730  0.47179    
price       -1.0444140  0.8343573  -1.252  0.22180    
temp         0.0034584  0.0004455   7.762  3.1e-08 ***
income       0.0033078  0.0011714   2.824  0.00899 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.001356648)

    Null deviance: 0.125523  on 29  degrees of freedom
Residual deviance: 0.035273  on 26  degrees of freedom
AIC: -107.24

Number of Fisher Scoring iterations: 2
m02 <- glm(cons~ price*temp + income, data = Icecream)
summary(m02)

Call:
glm(formula = cons ~ price * temp + income, data = Icecream)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.067043  -0.010332   0.001666   0.016298   0.073229  

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -0.054024   0.807061  -0.067    0.947  
price       -0.099346   2.976747  -0.033    0.974  
temp         0.008859   0.016312   0.543    0.592  
income       0.003209   0.001229   2.612    0.015 *
price:temp  -0.019675   0.059399  -0.331    0.743  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.001404748)

    Null deviance: 0.125523  on 29  degrees of freedom
Residual deviance: 0.035119  on 25  degrees of freedom
AIC: -105.37

Number of Fisher Scoring iterations: 2
library(texreg)
screenreg(list(m0, m01, m02))

================================================
                Model 1   Model 2      Model 3  
------------------------------------------------
(Intercept)       0.92 *     0.20        -0.05  
                 (0.40)     (0.27)       (0.81) 
price            -2.05      -1.04        -0.10  
                 (1.44)     (0.83)       (2.98) 
temp                         0.00 ***     0.01  
                            (0.00)       (0.02) 
income                       0.00 **      0.00 *
                            (0.00)       (0.00) 
price:temp                               -0.02  
                                         (0.06) 
------------------------------------------------
AIC             -75.25    -107.24      -105.37  
BIC             -71.05    -100.23       -96.96  
Log Likelihood   40.63      58.62        58.69  
Deviance          0.12       0.04         0.04  
Num. obs.        30         30           30     
================================================
*** p < 0.001, ** p < 0.01, * p < 0.05

From this we can see that model 1 shows no significance. model 2 has two significant predictors (temp and income). Model threeshows only income as a signifcant predictor. Model 2 is the best model.

anova(m0, m01, m02, test = "Chisq")
Analysis of Deviance Table

Model 1: cons ~ price
Model 2: cons ~ price + temp + income
Model 3: cons ~ price * temp + income
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1        28   0.117064                          
2        26   0.035273  2 0.081792 2.273e-13 ***
3        25   0.035119  1 0.000154    0.7405    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Again we see model two has the most significance.
library(visreg)
visreg(m02, "income", scale="response")
  Note that you are attempting to plot a 'main effect' in a model that contains an
  interaction.  This is potentially misleading; you may wish to consider using the 'by'
  argument.
Conditions used in construction of plot
price: 0.277
temp: 49.5

The higher the income the more consumption -makes sense.

visreg(m02, "income", by="temp", scale="response")

The higher the temperature, the more ice cream consumption especially when looking at income.

visreg(m02, "price", by="income", scale="response")

As the price increased, consumption decreased. Those who had a higher income started off with a higher consuption but it slowly decreased as well as price increased.

visreg(m02, "income", by="price", scale="response")

visreg(m0, "price", scale="response")

As the price increased, the ice cream consumption decreased.

visreg(m01, "temp", by="income", scale="response")

So I really wanted to see what would influence people to consume more ice cream. I really thought price would be much more signifcant but what I noticed is that the temperature and and a person’s income per week influenced how much someone consumed ice cream. Price did influence when compared to the income.

The more a person made the more ice cream they would eat but as the price increased, those who made more income were not consuming as much as I thought they would. It makes sense that as temperature rises, the more people want to be cooled and eat something refreshing. Something to note for future research : age, region / state and gender would be interesting to analyze. Also they looked at data between 03/18/1951 - 07/11/1953, I wonder why they didnt just start and end in January.

LS0tDQp0aXRsZTogIkhvbWV3b3JrIDUgTW9yZSBMb2cgUmVnIE1vZGVscyINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCiMjIyMjIEkgd2lsbCBiZSBhbmFseXppbmcgYSBkYXRhIGZyYW1lIGNvbnRhaW5pbmcgZGF0YSBvbiBUaGUgQ29uc3VtcHRpb24gb2YgSWNlIENyZWFtLiBUaGUgdmFyaWFibGVzIGluY2x1ZGU6DQoNCjEuY29uc3VtcHRpb24gb2YgaWNlIGNyZWFtIHBlciBoZWFkIGluIHBpbnRzDQoyLmluY29tZSBvZiB0aGUgYXZlcmFnZSBmYW1pbHkgcGVyIHdlZWsNCjMudGhlIHByaWNlIG9mIGljZSBjcmVhbSBwZXIgcGludA0KNC5hbmQgdGhlIGF2ZXJhZ2UgdGVtcGVyYXR1cmUgaW4gdGhlIFUuUy4NCg0KVGhlIG9ic2VydmF0aW9uIG51bWJlcnMgd2VyZSAzMCBhbmQgaXQgd2FzIGFjcm9zcyA0IHdlZWtzIC4NCg0KYGBge3J9DQpsaWJyYXJ5KEVjZGF0KQ0KbGlicmFyeShFY2Z1bikNCmBgYA0KDQoNCmBgYHtyfQ0KZGF0YSgiSWNlY3JlYW0iKQ0KaGVhZChJY2VjcmVhbSkNCmBgYA0KDQpgYGB7cn0NCmxpYnJhcnkoZ2dwbG90MikNCmdncGxvdChJY2VjcmVhbSwgYWVzKHg9Y29ucykpICsgZ2VvbV9iYXIoZmlsbD0iYmx1ZSIsIHdpZHRoPS4wMDQpDQpgYGANCiMjIyMjU28gaXQgc2VlbXMgYXMgdGhvdWdoIHRoZSBoaWdoZXN0IGNvbnN1bXB0aW9uIG9mIGljZSBjcmVhbSBpcyBhYm91dCAyIHBpbnRzLg0KDQoNCmBgYHtyfQ0KZ2dwbG90KEljZWNyZWFtLCBhZXMoeD1jb25zLCB5PXRlbXApKSArIGdlb21fYmFyKHN0YXQ9ImlkZW50aXR5IixmaWxsPSJibHVlIiwgd2lkdGg9LjAwNCkNCmBgYA0KIyMjIyBBbmQgdGhlbiB3ZSBjYW4gc2VlIGhlcmUgdGhhdCBhcyB0aGUgdGVtcGVyYXR1cmUgcmlzZXMsIHNvIGRvZXMgdGhlIGNvbnN1bXB0aW9uIG9mIGljZSBjcmVhbS4NCg0KDQpgYGB7cn0NCm0xPC1sbSAoY29uc34gcHJpY2UsIGRhdGEgPSBJY2VjcmVhbSkNCnN1bW1hcnkobTEpDQpgYGANCg0KYGBge3J9DQptMjwtbG0gKGNvbnMgfiBwcmljZSArIGluY29tZSwgZGF0YSA9IEljZWNyZWFtKQ0Kc3VtbWFyeShtMikNCmBgYA0KDQpgYGB7cn0NCm0zIDwtbG0gKGNvbnMgfiBwcmljZSArIHRlbXAsIGRhdGEgPSBJY2VjcmVhbSkNCnN1bW1hcnkobTMpDQpgYGANCiMjIyMgVGVtcGVyYXR1cmUgaGFzIGEgc3Ryb25nIHNpZ25pZmljYW5jZSBvbiBjb25zdW1wdGlvbiBpbiB0aGlzIGxpbmVhciBtb2RlbC4NCg0KDQoNCmBgYHtyfQ0KbGlicmFyeSh0ZXhyZWcpDQpzY3JlZW5yZWcobGlzdChtMSxtMixtMykpDQpgYGANCg0KYGBge3J9DQptNDwtbG0oY29ucyB+IHByaWNlKnRlbXAsIGRhdGE9SWNlY3JlYW0pDQpzdW1tYXJ5KG00KQ0KYGBgDQoNCg0KYGBge3J9DQptMCA8LSBnbG0oY29uc34gcHJpY2UsIGRhdGEgPSBJY2VjcmVhbSkNCnN1bW1hcnkobTApDQpgYGANCg0KDQpgYGB7cn0NCm0wMSA8LSBnbG0oY29uc34gcHJpY2UgKyB0ZW1wICsgaW5jb21lLCBkYXRhID0gSWNlY3JlYW0pDQpzdW1tYXJ5KG0wMSkNCmBgYA0KDQpgYGB7cn0NCm0wMiA8LSBnbG0oY29uc34gcHJpY2UqdGVtcCArIGluY29tZSwgZGF0YSA9IEljZWNyZWFtKQ0Kc3VtbWFyeShtMDIpDQpgYGANCg0KDQoNCmBgYHtyfQ0KbGlicmFyeSh0ZXhyZWcpDQpzY3JlZW5yZWcobGlzdChtMCwgbTAxLCBtMDIpKQ0KYGBgDQojIyMjIEZyb20gdGhpcyB3ZSBjYW4gc2VlIHRoYXQgbW9kZWwgMSBzaG93cyBubyBzaWduaWZpY2FuY2UuIG1vZGVsIDIgaGFzIHR3byBzaWduaWZpY2FudCBwcmVkaWN0b3JzICh0ZW1wIGFuZCBpbmNvbWUpLiBNb2RlbCB0aHJlZXNob3dzIG9ubHkgaW5jb21lIGFzIGEgc2lnbmlmY2FudCBwcmVkaWN0b3IuIE1vZGVsIDIgaXMgdGhlIGJlc3QgbW9kZWwuDQoNCg0KYGBge3J9DQphbm92YShtMCwgbTAxLCBtMDIsIHRlc3QgPSAiQ2hpc3EiKQ0KDQpgYGANCiMjIyMjIEFnYWluIHdlIHNlZSBtb2RlbCB0d28gaGFzIHRoZSBtb3N0IHNpZ25pZmljYW5jZS4NCmBgYHtyfQ0KbGlicmFyeSh2aXNyZWcpDQp2aXNyZWcobTAyLCAiaW5jb21lIiwgc2NhbGU9InJlc3BvbnNlIikNCg0KYGBgDQoNCiMjIyMgVGhlIGhpZ2hlciB0aGUgaW5jb21lIHRoZSBtb3JlIGNvbnN1bXB0aW9uIC1tYWtlcyBzZW5zZS4NCg0KYGBge3J9DQp2aXNyZWcobTAyLCAiaW5jb21lIiwgYnk9InRlbXAiLCBzY2FsZT0icmVzcG9uc2UiKQ0KYGBgDQojIyMjIFRoZSBoaWdoZXIgdGhlIHRlbXBlcmF0dXJlLCB0aGUgbW9yZSBpY2UgY3JlYW0gY29uc3VtcHRpb24gZXNwZWNpYWxseSB3aGVuIGxvb2tpbmcgYXQgaW5jb21lLiANCg0KYGBge3J9DQp2aXNyZWcobTAyLCAicHJpY2UiLCBieT0iaW5jb21lIiwgc2NhbGU9InJlc3BvbnNlIikNCmBgYA0KDQojIyMjIEFzIHRoZSBwcmljZSBpbmNyZWFzZWQsIGNvbnN1bXB0aW9uIGRlY3JlYXNlZC4gVGhvc2Ugd2hvIGhhZCBhIGhpZ2hlciBpbmNvbWUgc3RhcnRlZCBvZmYgd2l0aCBhIGhpZ2hlciBjb25zdXB0aW9uIGJ1dCBpdCBzbG93bHkgZGVjcmVhc2VkIGFzIHdlbGwgYXMgcHJpY2UgaW5jcmVhc2VkLg0KYGBge3J9DQoNCnZpc3JlZyhtMDIsICJpbmNvbWUiLCBieT0icHJpY2UiLCBzY2FsZT0icmVzcG9uc2UiKQ0KYGBgDQoNCg0KDQpgYGB7cn0NCnZpc3JlZyhtMCwgInByaWNlIiwgc2NhbGU9InJlc3BvbnNlIikNCmBgYA0KIyMjIyBBcyB0aGUgcHJpY2UgaW5jcmVhc2VkLCB0aGUgaWNlIGNyZWFtIGNvbnN1bXB0aW9uIGRlY3JlYXNlZC4NCmBgYHtyfQ0KdmlzcmVnKG0wMSwgInRlbXAiLCBieT0iaW5jb21lIiwgc2NhbGU9InJlc3BvbnNlIikNCmBgYA0KDQojIyMjIFNvIEkgcmVhbGx5IHdhbnRlZCB0byBzZWUgd2hhdCB3b3VsZCBpbmZsdWVuY2UgcGVvcGxlIHRvIGNvbnN1bWUgbW9yZSBpY2UgY3JlYW0uIEkgcmVhbGx5IHRob3VnaHQgcHJpY2Ugd291bGQgYmUgbXVjaCBtb3JlIHNpZ25pZmNhbnQgYnV0IHdoYXQgSSBub3RpY2VkIGlzIHRoYXQgdGhlIHRlbXBlcmF0dXJlIGFuZCBhbmQgYSBwZXJzb24ncyBpbmNvbWUgcGVyIHdlZWsgaW5mbHVlbmNlZCBob3cgbXVjaCBzb21lb25lIGNvbnN1bWVkIGljZSBjcmVhbS4gUHJpY2UgZGlkIGluZmx1ZW5jZSB3aGVuIGNvbXBhcmVkIHRvIHRoZSBpbmNvbWUuIA0KDQojIyMjVGhlIG1vcmUgYSBwZXJzb24gbWFkZSB0aGUgbW9yZSBpY2UgY3JlYW0gdGhleSB3b3VsZCBlYXQgYnV0IGFzIHRoZSBwcmljZSBpbmNyZWFzZWQsIHRob3NlIHdobyBtYWRlIG1vcmUgaW5jb21lIHdlcmUgbm90IGNvbnN1bWluZyBhcyBtdWNoIGFzIEkgdGhvdWdodCB0aGV5IHdvdWxkLiBJdCBtYWtlcyBzZW5zZSB0aGF0IGFzIHRlbXBlcmF0dXJlIHJpc2VzLCB0aGUgbW9yZSBwZW9wbGUgd2FudCB0byBiZSBjb29sZWQgYW5kIGVhdCBzb21ldGhpbmcgcmVmcmVzaGluZy4gU29tZXRoaW5nIHRvIG5vdGUgZm9yIGZ1dHVyZSByZXNlYXJjaCA6IGFnZSwgcmVnaW9uIC8gc3RhdGUgYW5kIGdlbmRlciB3b3VsZCBiZSBpbnRlcmVzdGluZyB0byBhbmFseXplLiBBbHNvIHRoZXkgbG9va2VkIGF0IGRhdGEgYmV0d2VlbiAwMy8xOC8xOTUxIC0gMDcvMTEvMTk1MywgSSB3b25kZXIgd2h5IHRoZXkgZGlkbnQganVzdCBzdGFydCBhbmQgZW5kIGluIEphbnVhcnkuIA0K