Abstract

My ME/CFS regularly improves during the summer. In a previous study, I found that the improvement correlates with air temperature but not with other environmental parameters. In a simulation of summerish conditions performed in a room, during cold months, by artificially increasing air temperature, infrared radiation, and air humidity, the correlation between my symptoms and air temperature tended to be negative. Here I report on a further study where I consider again environmental parameters from a few locations in the temperate and subtropical regions of both hemispheres and I search for the best regression model between them and my status (n=20). The best regression model, among the ones tested, is Y=\beta_0+\beta_1 t +\beta_2 T^4 where Y is my status (the higher the value, the less the symptoms), t is air temperature (°C) and T is absolute air temperature (K). This regression model approximates the flux of thermic energy between my body and the environment, when indoors, but with two different signs for the coefficients of t and of T^4. I have used then the best regression model to predict the score over the year in several cities. These predictions can be used to find the optimal spot to live in, as a function of the month of the year.

Introduction

My ME/CFS improves during summer, in the period of the year that goes from May/June to the end of September (when I am in Rome, Italy). I don’t know why. In a previous study (Maccallini P. 2022), I analyzed the correlation between my symptoms and several environmental parameters, during three periods covering a total of 190 days: in one, I improved while the season was going toward summer, in another one we were in full-blown summer, and in the third one my functioning declined while summer was fading away. The statistical analysis showed a positive correlation between my functionality score and the temperature of the air and a negative correlation between the score and air density and dry air density (remember that the higher the functionality score, the better I feel). No correlation was found between my functionality score and: relative humidity, absolute humidity, atmospheric pressure, and the concentration of particulate with a diameter below or equal to 2.5 \mu m. In another study, I simulated summerish conditions in a room for five months, during cold months, by means of a heat pump, an infrared lamp, an air moisturizer, and mylar panels (Maccallini P. 2022) and despite an improvement during the experiment, the correlation between my status and air temperature tended to be negative. In the present study, I consider again environmental parameters from a few locations in the temperate and subtropical regions of both hemispheres and I search for the best regression model between them and my status.

Figure 1. The duplets City-Month are reported in the euclidean space Cl, R, T. The points City-Months that have already been classified based on the score I registered for them, are colored green (score of 6), red (score of 4), yellow with a black contour line (score of 5), and yellow with a green contour line (score of 5.5). Black dots are doublets City-Month without a classification. CO = Corrientes, Argentina; AS = Asunción, Paraguay; SM = Santa Maria, Cape Verde. F = February; Ma = March, and so on.

Methods and Results

I have considered a set of 20 points defined by the two coordinates City and Month. To each of them, I have assigned a score, from 4 to 6, based on the degree of severity (the lower the score, the worse the disease) of my symptoms in that particular city, for that specific month. For each point City-Month I have collected the monthly mean air temperature (T), mean cloudiness (Cl), and mean solar radiation (R) (Figure 1). I have considered several models of regression between the score (response variable) and the three explanatory variables T, Cl, and R (Table 1). For each regression model, I have calculated the p values of the estimates of the coefficients, and I have also calculated the p value of the F test, given by

(1) \,\,\,p\,=\,P\left(TS>\frac{R^2 \left(n-r-1\right)}{\left(1-R^2\right)r}\right)

where TS follows an F law of parameters r and n-r-1, with r the number of explanatory variables and n the number of data points. All the calculations and the plots of this article (except for Figure 3) have been performed by two custom scripts (reported at the end), one in Octave, and one in R. For the script in Octave I wrote a subroutine for the F-test. The R script also calculates Spearman’s, and Pearson’s correlation coefficients between the variables, including their p values and the corrections for temperature (Table 2).

Model\beta_0\beta_1\beta_2\beta_3F test
1)\,Y=\beta_0+\beta_1 R +\beta_2 t1.90**0.1100.129*–1.98e-4
2)\,Y=\beta_0+\beta_1 Cl +\beta_2 t0.5140.0210.198**–2.14e-4
3)\,Y=\beta_0+\beta_1 R +\beta_2 Cl + \beta_3 t0.2080.1360.02700.164*5.78e-4
4)\,Y=\beta_0+\beta_1 R +\beta_2 T^4-5.540.1001.36e-9*–1.51e-4
5)\,Y=\beta_0+\beta_1 R +\beta_2 t + \beta_3 T^4-1.49e2*0.131-2.60*2.77e-8*3.58e-5
6)\,Y=\beta_0+\beta_1 T^4-7.38**1.70e-9***––2.85e-5
7)\,Y=\beta_0+\beta_1 t +\beta_2 T^4-1.44e2*-2.47*2.69e-8*–1.28e-5

Table 1. For each regression model, I have indicated the values of the estimates of the coefficients (with their degree of significance), and the p value of the F test. The meaning of the asterisks is * p < 0.05; ** p<0.001; *** p = 0. T = mean daily absolute temperature; t = meant daily temperature in Celsius, R = mean daily solar radiation; Cl = fraction of the day mainly cloudy.

I have used then the best regression model to predict the score along the year in several cities in the temperate and subtropical regions of both hemispheres (Figure 2). The diagrams have been interpolated by using cubic splines, to make them more smooth. These predictions can be used to find the best spot to live in, as a function of the month of the year.

Figure 2. Prediction for the score by the regression model Y=\beta_0+\beta_1 T^4. RM = Rome; AR = Arona, Tenerife, Spain; RO = Rosario, Argentina; BA = Buenos Aires, Argentina; AS = Asunción, Paraguay; CO = Corrientes, Argentina; SM = Santa Maria, Cape Verde.

Cor. between
Y and		Spearman
(p value)		Pearson
(p value)		Pearson
corrected for t
(p value)		Pearson
corrected for R
(p value)
R0.758
(1.10e-4)		0.706
(0.509e-4)		0.204
(0.401)		–
Cl-0.511
(0.214e-2)		-0.583
(0.696e-3)		0.182
(0.456)		–
t
0.875
(4.49e-7)		0.786
(4.04e-5)		–0.519
(2.27e-2)
T^40.875
(4.49e-7)		0.794
(2.85e-5)		0.515
(2.43e-2)		–

Table 2. Correlation coefficients (Spearman and Person) between the variables. Corrected Pearson correlations were also computed, when possible. T = mean daily absolute temperature; t = meant daily temperature in Celsius, R = mean daily solar radiation; Cl = fraction of the day mainly cloudy.

I also used the 198 data points from my previous study (Maccallini P. 2022) to test models 7 and 6 of regression. I considered the mean data on three days for all the variables (Table 3). In this case, model 6 performs better. We still have a good level of significance in the F test, and we still have a negative coefficient for t in model 7. If we use this model to make previsions on the score as a function of the month in various cities, we get Figure 3. Consider that this regression model is based on a scoring system that is slightly lower if compared with the one used for the previous analysis.

Model\beta_0\beta_1\beta_2\beta_3F test
6)\,Y=\beta_0+\beta_1 T^48.19e-15.31e-10***––1.54e-6
7)\,Y=\beta_0+\beta_1 t +\beta_2 T^4-4.67e1.-9.09e-1.9.41e-9.–2.14e-6

Table 3. For each regression model, I have indicated the values of the estimates of the coefficients (with their degree of significance), and the p value of the F test. The meaning of the asterisks is: . p<01, * p < 0.05; ** p<0.001; *** p = 0. T = mean daily absolute temperature; t = meant daily temperature in Celsius, R = mean daily solar radiation; Cl = fraction of the day mainly cloudy.

Figure 2. Prediction for the score by the regression model 6 of Table 3. RM = Rome; AR = Arona, Tenerife, Spain; RO = Rosario, Argentina; BA = Buenos Aires, Argentina. AS = Asunción, Paraguay; CO = Corrientes, Argentina; SM = Santa Maria, Cape Verde.

Discussion

The best regression model is Y=\beta_0+\beta_1 t +\beta_2 T^4, followed by the model \beta_0+\beta_1 T^4, and model Y=\beta_0+\beta_1 R +\beta_2 t + \beta_3 T^4. These three models have a physical meaning; namely, they express the flux of thermic energy between my body and the environment. The human body must balance the energy exchanged with the environment to maintain thermic homeostasis. The power P exchanged with the environment can be written as

(2) \,\,\,P\,=\,R+CV+CD+E

where R is the power exchange in the form of radiant energy, C is the exchange by conduction, CV is the exchange by convection, and E is the dispersion of power by evaporation. This is a simplification of complex phenomena that should also take into account the exchange of thermic power associated with the passage of atmospheric gas through the respiratory system. The power R can be further specified in the form R\,=\,K_1 R_{sun} + K_2(T_me^4-T_env^4) where K_1 R_{sun} is the power absorbed from the Sun, and the other addend is the radiant power exchanged with the environment, whit K_2 a constant that expresses radiant and geometrical properties of both my body and of the environment. While R_{sun} is radiant power in the form of short wave radiation (wavelength below 3 \mu), the other addend represents radiant power exchanged by bodies at low temperatures with a wavelength above 3 \mu (according to the Stephan-Boltzman law) (R). But what temperature is T_{env}? One first approximation might be the absolute temperature of the air, because in warm months, with the windows open, the objects surrounding us tend to acquire a superficial temperature similar to that of the air. To show that point, I analyzed the power emitted by the soil in Rome, during the year 2020, and the temperature of the air (data from ARPALAZIO). If we indicate J the daily mean power emitted per surface unit by the soil in the hemispace, the Stefan-Boltzmann law states that

(3) \,\,\,J\,=\,a\chi T_s^4

where a is a constant below 1 dependent on the material of the surface of the soil, \chi\,=\,5.6696\cdot10^{-8}\frac{W}{m^2 K^4} is the Stefan-Boltzmann constant, and T_s is the absolute temperature of the soil. If the assertion that the temperature of the soil is about the same as the temperature of the air was correct, the regression for the following linear model should have a very low p value of the F test:

(4) \,\,\,J\,=\,B_0 + B_1(273.15 + t_a)^4

And this is in fact the case, we have the regression in Figure 3, with a p value of the F test below 2.2\cdot 10^{-16}. For further details on these calculations, see (Maccallini P, 2022).

Figure 3. Linear regression between the mean daily irradiance from the Earth and the fourth power of the absolute temperature of the air (daily mean value). Data from ARPA, Rome, 2020. The analysis is described in detail in (Maccallini P, 2022).

As for the exchange of power by conduction and convection, in both cases, they are expressed by a constant multiplied by the difference in temperature between my skin and the air. That said, we can write again eq. 1 as follows:

(5) \,\,\,P\,=\,K_1 R_{sun}+K_2 (T_{skin}^4-T_{air}^4)+K_3 (t_{skin} - t_{air})

were we have removed the term linked to evaporation. In other words, the power exchanged by my body with the environment could be expressed, in first approximation, by the equation

(6) \,\,\,P\,=\,B_0+B_1 R_{sun}+B_2 t_{air} + B_3 T_{air}^4

This is model 5 of Table 1. If we consider the case of a subject who rarely leaves home, we can get rid of R_{sun} and we get model 7, the best model of regression. In a previous article, I showed that the best correlation between my status and several environmental parameters was with air temperature (R), and yet in an experiment where I simulated summer in a room with the use of a heat pump, an infrared lamp, and an air moisturizer, there was no correlation between my status and air temperature (R) (in fact the correlation had a tendency to be negative). From the perspective of this new study, I would explain the previous results as follows: in the first study performed during worm months using environmental parameters we caught the Spearman’s correlation between my status and air temperature which was a reflection of the correlation between my status and the longwave infrared radiation from the environment (proportional to the fourth power of the absolute temperature of the air); in the second study, performed during cold months within a room with artificially heated air, there was no correlation between air temperature and the infrared emission of the environment (which came from a lamp), so we found a correlation that tended to be negative perhaps because when I was feeling worse I raised the temperature of the air. It should be noted though that in models 7 and 5 the coefficient of t is negative, therefore it seems that while the longwave infrared radiation plays a positive role, the temperature of the air does not.

Conclusion

The three best regression models suggest that my improvement is linked to the thermal energy balance of the body: the lower the dispersion of radiant thermic energy, the better I feel. But at the same time, the temperature of the air has a detrimental effect. But other interpretations might be possible: there might be an effect of thermal radiation on my biology that is not linked to thermic homeostasis but instead to some other biological targets. Moreover, there might be other models of regression that I have not considered, that perform better than the ones mentioned.

Scripts

The first script is written in Octave. It calculates the regressions and their F tests (by using a custom subroutine) and it plots Figure 1 and Figure 2. It also writes the data in csv format for the second script (in R), which performs the same statistical analyses and also calculates the p values associated with each coefficient estimate.

% file name = Spot_finder_2
% date of creation = 29/01/2023
%
close all
clear all
#
pkg load statistics % we need this package for the F law
spline = 1 % it is one if we want a cubic spline interpolation for the last plot
#
#-------------------------------------------------------------------------------
# Data
#-------------------------------------------------------------------------------
#
months = ["Je"; "F"; "Mr"; "Ap"; "Ma"; "Jn"; "Jl"; "Au"; "S"; "O"; "N"; "D"];
%
% data of Rome (RM), Italy
%
k=1;
Spot(1,:,k) = [2,3,4.3,5.6,6.8,7.5,7.5,6.5,4.9,3.4,2.2,1.8]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [47,44,44,43,39,25,13,18,32,43,47,46]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [7,8,11,13,18,22,25,25,21,17,12,8]; % Mean air temperature (°C)
Spot(4,:,k) = [3,3,6,8,12,16,18,18,15,12,7,4]; % Min air temperature (°C)
Spot(5,:,k) = [12,13,16,19,23,27,31,31,27,22,17,13]; % Max air temperature (°C)
%
% data of Arona (AR), Tenerife, Spain
%
k=2;
Spot(1,:,k) = [3.9,4.9,6.1,7.2,7.8,8.2,8.0,7.4,6.3,5.1,4.0,3.5]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [32,28,26,24,21,9,2,6,21,34,38,36]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [19,19,19,20,21,22,24,25,25,24,22,20]; % Mean air temperature (°C)
Spot(4,:,k) = [16,16,16,17,18,20,21,22,22,21,19,17]; % Min air temperature (°C)
Spot(5,:,k) = [22,22,23,23,24,26,28,28,28,27,25,23]; % Max air temperature (°C)
%
% data of Rosario (RO), Santa Fe, Argentina
%
k=3;
Spot(1,:,k) = [7.6,6.7,5.6,4.3,3.2,2.7,3.0,3.9,5.1,6.3,7.3,7.7]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [28,28,28,34,45,47,44,40,36,32,30,28]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [25,23,22,18,14,11,10,12,15,18,21,24]; % Mean air temperature (°C)
Spot(4,:,k) = [19,18,17,13,10,7,6,7,9,13,16,18]; % Min air temperature (°C)
Spot(5,:,k) = [30,29,27,23,19,16,16,18,21,24,27,29]; % Max air temperature (°C)
%
% data of Buenos Aires (BA), Argentina
%
k=4;
Spot(1,:,k) = [7.6,6.7,5.5,4.1,3.0,2.5,2.7,3.6,4.8,6.2,7.3,7.8]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [29,31,31,38,48,50,48,45,41,37,33,30]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [24,23,22,18,15,12,11,13,14,17,20,23]; % Mean air temperature (°C)
Spot(4,:,k) = [21,20,19,15,12,9,9,10,11,14,17,19]; % Min air temperature (°C)
Spot(5,:,k) = [28,27,25,21,18,15,14,16,18,21,24,27]; % Max air temperature (°C)
%
% data of Asuncion (AS), Paraguay
%
k=5;
Spot(1,:,k) = [7.0,6.6,6.0,5.0,4.1,3.5,3.8,4.5,5.3,6.1,6.9,7.1]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [45,42,33,31,32,37,35,32,32,38,38,42]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [28,27,26,23,20,18,18,20,21,24,25,27]; % Mean air temperature (°C)
Spot(4,:,k) = [23,23,21,19,16,14,14,15,16,19,21,22]; % Min air temperature (°C)
Spot(5,:,k) = [33,32,31,28,25,23,23,25,27,29,31,32]; % Max air temperature (°C)
%
% data of Corrientes (CO), Argentina
%
k=6;
Spot(1,:,k) = [7.2,6.7,5.9,4.8,3.9,3.4,3.6,4.3,5.3,6.2,7.0,7.3]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [39,35,30,30,33,37,35,33,32,34,34,36,]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [27,26,25,22,18,16,15,17,19,22,24,26]; % Mean air temperature (°C)
Spot(4,:,k) = [22,22,20,17,14,12,11,12,14,17,19,21]; % Min air temperature (°C)
Spot(5,:,k) = [33,32,30,26,23,21,21,23,25,27,29,31]; % Max air temperature (°C)
%
% data of SantaMaria (SM), Capo Verde
%
k=7;
Spot(1,:,k) = [5.0,5.9,6.8,7.3,7.3,7.2,6.7,6.2,5.8,5.6,5.0,4.6]; % Mean daily incident solar radiation (kWh/m^2)
Spot(2,:,k) = [42,35,34,33,25,19,36,53,56,48,54,49]; % Fraction of the day mainly cloudy (%)
Spot(3,:,k) = [22,22,22,22,23,24,25,27,27,27,25,23]; % Mean air temperature (°C)
Spot(4,:,k) = [20,19,20,20,21,22,23,25,25,24,23,21]; % Min air temperature (°C)
Spot(5,:,k) = [25,25,25,26,26,27,28,30,30,30,28,26]; % Max air temperature (°C)
%
names = {"RM","AR","RO","BA","AS","CO","SM"};
types = {"-k","--k",":k","-r","--r",":r","-b","--b","-.b",}
%
for i=1:5
  max_measure(i) = max([Spot(i,:,1),Spot(i,:,2),Spot(i,:,3),Spot(i,:,4),Spot(i,:,5),Spot(i,:,6),Spot(i,:,7)])
endfor
%
%-------------------------------------------------------------------------------
% Subroutine for F-test
%-------------------------------------------------------------------------------
%
function [pval,R2,TS] = F_test (M,b,beta)
  r = length(M(1,:))-1; % number of explanatory variables
  n = length(M(:,1)) % number of data points
  R = b-M*beta; % residues
  SSR = 0;
  Sbb = 0;
  for i=1:n
    SSR = SSR + R(i)^2; % sum of squared residues
    Sbb = Sbb + b(i)^2;
  endfor
  Sbb = Sbb - n*(mean(b))^2;
  R2 = (Sbb - SSR)/Sbb; % R squared
  TS = R2*(n-r-1)/(r*(1-R2)); % test statistics
  pval = 1 - fcdf(TS,r,n-r-1); % F-test
endfunction
%
%-------------------------------------------------------------------------------
% Subroutine plot_score
%-------------------------------------------------------------------------------
%
function Plot_score (score, h, model, names, types, spline)
    figure(5) % Plotting of the Scoring for the query spots
    %
    for k=1:length(names)
      if (spline==1)
        yi = interp1([1:1:12],score(:,k),[1:0.25:12],"spline")
        plot([1:0.25:12],yi,char(types{k}),"linewidth",1.5)
        hold on
      else
        plot([1:1:12],score(:,k),char(types{k}),"linewidth",1.5)
        hold on
      endif
    endfor
    plot([1,12],[5,5],"--g","linewidth",1.5)
    title(model{h})
    legend(char(names{1:length(names)}),"location","southwest",'fontsize',15)
    xlabel("Month",'fontsize',15);
    ylabel("Score",'fontsize',15);
    grid on
    grid minor on
endfunction
%
%-------------------------------------------------------------------------------
figure 1 % mean temp, irradiation, cloudiness with storing of data for regressions
%-------------------------------------------------------------------------------
%
% plotting Rome
%
k=1;
n=1;
for i=1:12
  if (and(i>=6,i<9))
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
    b(n)=6;
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  if (or(i<6,i>9))
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','r','markeredgecolor','k')
  endif
  if (or(and(i>3,i<5),i>9))
    b(n)=4;
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  if (i==5)
    b(n)=5;
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  if (i==9)
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','y','markeredgecolor','g')
    b(n)=5.7;
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  hold on
endfor
%
% plotting Arona
%
k=2;
for i=1:12
  if (i==5)
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
    b(n)=6
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  if (and(i>=1,i<5))
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','r','markeredgecolor','k')
    b(n)=4
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
    hold on
endfor
%
% plotting Rosario
%
k=3;
for i=1:12
  if (and(i>=1,i<3))
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
    b(n)=6
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  if (i==3)
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','y','markeredgecolor','k')
    b(n)=5.5
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  hold on
endfor
%
% plotting CABA
%
k=4;
for i=1:12
  if (or(i==12,i==1,i==2))
    plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
    b(n)=6
    M(n,1:4)=[1,Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)];
    n=n+1
  endif
  if (i<12)
    hold on
  endif
endfor
%
% plotting query spots
%
% Asuncion
%
k=5;
for i=2:4
  plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
% Corrientes
%
k=6;
for i=2:4
  plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
% Santa Maria
%
k=7;
for i=2:5
  plot3(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(1,i,k),Spot(2,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
grid on
grid minor on
xlabel("Mean daily incident solar radiation (kWh/m^{2})",'fontsize',15);
ylabel("Fraction of the day mainly cloudy (%)",'fontsize',15);
zlabel("Mean air temperature (°C)",'fontsize',15);
axis ('square')
axis([4,max_measure(1),0,max_measure(2),15,max_measure(3)])
%
%-------------------------------------------------------------------------------
% Regressions
%-------------------------------------------------------------------------------
%
n=n-1; % number of measures
%
b = transpose(b)
x = cell() % it contains the coefficients of the models
model = cell() % it contains the description of the models
matrix = cell() % it contains the matrix of the regression
k = 0;
%
for i=1:n
  A(i,1:3)=[1,M(i,2),M(i,4)]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*Rad + b2*t"
matrix{k} = A;
clear A
%
for i=1:n
  A(i,1:3)=[1,M(i,3),M(i,4)]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*Cl + b2*t"
matrix{k} = A;
clear A
%
for i=1:n
  A(i,1:4)=[1,M(i,2),M(i,3),M(i,4)]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*Rad + b1*Cl + b2*t"
matrix{k} = A;
clear A
%
for i=1:n
  A(i,1:3)=[1,M(i,2),(273.15+M(i,4))^4]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*Rad + b2*T^{4}"
matrix{k} = A;
clear A
%
for i=1:n
  A(i,1:4)=[1,M(i,2),M(i,4),(273.15+M(i,4))^4]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*Rad + b2*t + b3*T^{4}"
matrix{k} = A;
clear A
%
for i=1:n
  A(i,1:2)=[1,(273.15+M(i,4))^4]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*T^{4}"
matrix{k} = A;
clear A
%
for i=1:n
  A(i,1:3)=[1,M(i,4),(273.15+M(i,4))^4]
endfor
k = k+1;
[x{k}, std, mse, S] = lscov (A, b)
[pval(k), R2(k), TS(k)] = F_test (A, b, x{k})
model{k} = "socre = b0 + b1*t + b2*T^{4}"
matrix{k} = A;
clear A
%
%-------------------------------------------------------------------------------
% Selection of the best model
%-------------------------------------------------------------------------------
%
for i=1:length(pval)
  if (pval(i)==min(pval))
    best=i;
  endif
endfor
%
%-------------------------------------------------------------------------------
% Region of the plane Rad-Temp with high score according to the regression
% model
%-------------------------------------------------------------------------------
%
R(1)= 4;
R(2)= max_measure(1);
Temp_R(1) = ( 5.5 - x{1}(1) - x{1}(2)*R(1) )/x{1}(3);
Temp_R(2) = ( 5.5 - x{1}(1) - x{1}(2)*R(2) )/x{1}(3);
%
%-------------------------------------------------------------------------------
% Region of the plane Cloudiness-Temp with high score according to the regression
% model
%-------------------------------------------------------------------------------
%
Cl(1)= 0;
Cl(2)= max_measure(2);
Temp_Cl(1) = ( 5.5 - x{2}(1) - x{2}(2)*Cl(1) )/x{2}(3);
Temp_Cl(2) = ( 5.5 - x{2}(1) - x{2}(2)*Cl(2) )/x{2}(3);
%
%-------------------------------------------------------------------------------
figure 2 % mean temp and irradiation
%-------------------------------------------------------------------------------
%
% plotting Rome
%
k=1;
for i=1:12
  if (and(i>=6,i<9))
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (or(i<6,i>9))
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','r','markeredgecolor','k')
  endif
  if (i==9)
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','y','markeredgecolor','g')
  endif
  hold on
endfor
%
% plotting Arona
%
k=2;
for i=1:12
  if (i==5)
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (and(i>=1,i<5))
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','r','markeredgecolor','k')
  endif
    hold on
endfor
%
% plotting Rosario
%
k=3;
for i=1:12
  if (and(i>=1,i<3))
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (i==3)
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','y','markeredgecolor','k')
  endif
  hold on
endfor
%
% plotting CABA
%
k=4;
for i=1:12
  if (or(i==12,i==1,i==2))
    plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (i<12)
    hold on
  endif
endfor
%
% plotting query spots
%
% Asuncion
%
k=5;
for i=2:4
  plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(1,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
% Corrientes
%
k=6;
for i=2:4
  plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(1,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
% Santa Maria
%
k=7;
for i=3:5
  plot(Spot(1,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(1,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
plot([R(1),R(2)],[Temp_R(1),Temp_R(2)],"--k","linewidth",1)
%
grid on
grid minor on
xlabel("Mean daily incident solar radiation (kWh/m^{2})",'fontsize',15);
ylabel("Mean air temperature (°C)",'fontsize',15);
axis([R(1),R(2),15,max_measure(3)])
%
%-------------------------------------------------------------------------------
figure 3 % mean temp and cloudiness
%-------------------------------------------------------------------------------
%
% plotting Rome
%
k=1;
for i=1:12
  if (and(i>=6,i<9))
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (or(i<6,i>9))
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','r','markeredgecolor','k')
  endif
  if (i==9)
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','y','markeredgecolor','g')
  endif
  hold on
endfor
%
% plotting Arona
%
k=2;
for i=1:12
  if (i==5)
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (and(i>=1,i<5))
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','r','markeredgecolor','k')
  endif
    hold on
endfor
%
% plotting Rosario
%
k=3;
for i=1:12
  if (and(i>=1,i<3))
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (i==3)
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','y','markeredgecolor','k')
  endif
  hold on
endfor
%
% plotting CABA
%
k=4;
for i=1:12
  if (or(i==12,i==1,i==2))
    plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','g','markeredgecolor','k')
  endif
  if (i<12)
    hold on
  endif
endfor
%
% plotting query spots
%
% Buenos Aires
%
i=2;
k=4;
plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
text(Spot(2,i,k),Spot(3,i,k)-0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
%
% Asuncion
%
k=5;
for i=2:4
  plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(2,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
% Corrientes
%
k=6;
for i=2:4
  plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(2,i,k),Spot(3,i,k)-0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
% Santa Maria
%
k=7;
for i=3:5
  plot(Spot(2,i,k),Spot(3,i,k),'o','markersize',6,'markerfacecolor','k','markeredgecolor','k')
  text(Spot(2,i,k),Spot(3,i,k)+0.5, strcat(char(names{k}),", ",months(i,:)),'fontsize',15)
endfor
%
plot([Cl(1),Cl(2)],[Temp_Cl(1),Temp_Cl(2)],"--k","linewidth",1)
%
grid on
grid minor on
xlabel("Fraction of the day mainly cloudy (%)",'fontsize',15);
ylabel("Mean air temperature (°C)",'fontsize',15);
axis([Cl(1),Cl(2),15,max_measure(3)])
%
%-------------------------------------------------------------------------------
% Plotting the scoring for each city according to the best model
%-------------------------------------------------------------------------------
%
switch best
  case 1
    %
    % Scoring for the query spots according to model h
    %
    h=1;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*Spot(1,i,k)+x{h}(3)*Spot(3,i,k);
      endfor
    endfor
    %
    figure(5) % Plotting of the Scoring for the query spots
    %
    for k=1:length(names)
      plot([1:1:12],score(:,k),char(types{k}),"linewidth",1.5)
      hold on
    endfor
    plot([1,12],[5,5],"--g","linewidth",1.5)
    title(model{h},'fontsize',15)
    legend(char(names{1:length(names)}),"location","southwest",'fontsize',15)
    xlabel("Month",'fontsize',15);
    ylabel("Score",'fontsize',15);
    grid on
    grid minor on
  case 2
    %
    % Scoring for the query spots according to model h 
    %
    h=2;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*Spot(2,i,k)+x{h}(3)*Spot(3,i,k);
      endfor
    endfor
    %
    figure(5) % Plotting of the Scoring for the query spots
    %
    for k=1:length(names)
      plot([1:1:12],score(:,k),char(types{k}),"linewidth",1.5)
      hold on
    endfor
    plot([1,12],[5,5],"--g","linewidth",1.5)
    title(model{h},'fontsize',15)
    legend(char(names{1:length(names)}),"location","southwest",'fontsize',15)
    xlabel("Month",'fontsize',15);
    ylabel("Score",'fontsize',15);
    grid on
    grid minor on
  case 3
    %
    % Scoring for the query spots according to model h
    %
    h=3;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*Spot(1,i,k)+x{h}(3)*Spot(2,i,k)+x{h}(4)*Spot(3,i,k);
      endfor
    endfor
    %
    figure(5) % Plotting of the Scoring for the query spots
    %
    for k=1:length(names)
      plot([1:1:12],score(:,k),char(types{k}),"linewidth",1.5)
      hold on
    endfor
    plot([1,12],[5,5],"--g","linewidth",1.5)
    title(model{h},'fontsize',15)
    legend(char(names{1:length(names)}),"location","southwest",'fontsize',15)
    xlabel("Month",'fontsize',15);
    ylabel("Score",'fontsize',15);
    grid on
    grid minor on
  case 4
    %
    % Scoring for the query spots according to model h
    %
    h=4;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*Spot(1,i,k)+x{h}(3)*(Spot(3,i,k)^4);
      endfor
    endfor
    %
    figure(5) % Plotting of the Scoring for the query spots
    %
    for k=1:length(names)
      plot([1:1:12],score(:,k),char(types{k}),"linewidth",1.5)
      hold on
    endfor
    plot([1,12],[5,5],"--g","linewidth",1.5)
    title(model{h},'fontsize',15)
    legend(char(names{1:length(names)}),"location","southwest",'fontsize',15)
    xlabel("Month",'fontsize',15);
    ylabel("Score",'fontsize',15);
    grid on
    grid minor on
  case 5
    %
    % Scoring for the query spots according to model h
    %
    h=5;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*Spot(1,i,k)+x{h}(3)*Spot(3,i,k)+x{h}(4)*(Spot(3,i,k)^4);
      endfor
    endfor
    %
    figure(5) % Plotting of the Scoring for the query spots
    %
    for k=1:length(names)
      plot([1:1:12],score(:,k),char(types{k}),"linewidth",1.5)
      hold on
    endfor
    plot([1,12],[5,5],"--g","linewidth",1.5)
    title(model{h},'fontsize',15)
    legend(char(names{1:length(names)}),"location","southwest",'fontsize',15)
    xlabel("Month",'fontsize',15);
    ylabel("Score",'fontsize',15);
    grid on
    grid minor on
  case 6
    %
    % Scoring for the query spots according to model h
    %
    h=6;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*((273.15+Spot(3,i,k))^4);
      endfor
    endfor
    %
    Plot_score (score, h, model, names, types, spline)
  case 7
    %
    % Scoring for the query spots according to model h
    %
    h=7;
    for i=1:12
      for k=1:length(names)
        score(i,k) = x{h}(1)+x{h}(2)*Spot(3,i,k)+x{h}(3)*((273.15+Spot(3,i,k))^4);
      endfor
    endfor
    %
    Plot_score (score, h, model, names, types, spline)
endswitch
%
%-------------------------------------------------------------------------------
% We save the data for the regressions
%-------------------------------------------------------------------------------
%
save("M.csv","M")
save("Y.csv","b")
# file name: spot_finder
#
# Asunciòn, 9 marzo 2023
#
#-------------------------------------------------------------------------------
# This script performs the linear regressions performed also by the Octave script
# that goes under the same name
#-------------------------------------------------------------------------------
#
library(ppcor) # we need this for the partial correlation analysis
#
#-------------------------------------------------------------------------------
# Data 
#-------------------------------------------------------------------------------
#
Y<-read.csv(file = "Y.csv", header = F, sep =" ", comment.char = "#") # the response variable
Y<-Y[,2]
M<-read.csv(file = "M.csv", header = F, sep =" ", comment.char = "#")
M<-M[,2:5]
colnames(M)<-c("unity","Rad","Cl","meanT")
AT<-273.15+M$meanT # absolute temperature
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*R + x3*t) 
#-------------------------------------------------------------------------------
#
R_t_re<-lm(Y~M$Rad+M$meanT)
summary(R_t_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*Cl + x3*t) 
#-------------------------------------------------------------------------------
#
C_t_re<-lm(Y~M$Cl+M$meanT)
summary(C_t_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*R + x3*Cl + x4*t) 
#-------------------------------------------------------------------------------
#
R_C_t_re<-lm(Y~M$Rad+M$Cl+M$meanT)
summary(R_C_t_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*R + x4*T^4) 
#-------------------------------------------------------------------------------
#

R_T4_re<-lm(Y~M$Rad+I(AT^4))
summary(R_T4_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*R + x3*t + x4*T^4) 
#-------------------------------------------------------------------------------
#
R_t_T4_re<-lm(Y~M$Rad+M$meanT+I(AT^4))
summary(R_t_T4_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*T^4) 
#-------------------------------------------------------------------------------
#
T4_re<-lm(Y~I(AT^4))
summary(T4_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*t + x3T^4) 
#-------------------------------------------------------------------------------
#
t_T4_re<-lm(Y~M$meanT+I(AT^4))
summary(t_T4_re)
#
#-------------------------------------------------------------------------------
# Spearman coefficients and corresponding p values
#-------------------------------------------------------------------------------
#
Spearman_Rad<-cor.test(Y,M$Rad,alternative="two.sided",method="spearman",conf.level=0.95)
Spearman_Rad
Spearman_Cl<-cor.test(Y,M$Cl,alternative="two.sided",method="spearman",conf.level=0.95)
Spearman_Cl
Spearman_meanT<-cor.test(Y,M$meanT,alternative="two.sided",method="spearman",conf.level=0.95)
Spearman_meanT
#
#-------------------------------------------------------------------------------
# Pearson coefficients, corresponding p values, and partial correlations
#-------------------------------------------------------------------------------
#
Ro_Y_R<-cor.test(Y,M$Rad,alternative="two.sided",method="pearson",conf.level=0.95)
Ro_Y_R
Ro_Y_C<-cor.test(Y,M$Cl,alternative="two.sided",method="pearson",conf.level=0.95)
Ro_Y_C
Ro_Y_t<-cor.test(Y,M$meanT,alternative="two.sided",method="pearson",conf.level=0.95)
Ro_Y_t
Ro_Y_T4<-cor.test(Y,I(AT^4),alternative="two.sided",method="pearson",conf.level=0.95)
Ro_Y_T4
#
Ro_Y_R_t<-pcor.test(Y,M$Rad,M$meanT,method="pearson")
Ro_Y_R_t
Ro_Y_R_T4<-pcor.test(Y,M$Rad,I(AT^4),method="pearson")
Ro_Y_R_T4
Ro_Y_C_t<-pcor.test(Y,M$Cl,M$meanT,method="pearson")
Ro_Y_C_t
Ro_Y_C_T4<-pcor.test(Y,M$Cl,I(AT^4),method="pearson")
Ro_Y_C_T4
Ro_Y_t_R<-pcor.test(Y,M$meanT,M$Rad,method="pearson")
Ro_Y_t_R
Ro_Y_T4_R<-pcor.test(Y,I(AT^4),M$Rad,method="pearson")
Ro_Y_T4_R
Ro_Y_T4_T<-pcor.test(Y,I(AT^4),M$meanT,method="pearson")
Ro_Y_T4_T
#
#-------------------------------------------------------------------------------
# Data of the first environmental study
#-------------------------------------------------------------------------------
#
mydata<-read.csv(file = "correlation.csv", header = T, sep =";")
#
attach(mydata)
n<-length(Station)
Y<-c()
for (i in 1:n) Y[i]<-mean(c(Score1[i], Score2[i], Score3[i]), na.rm = T )
meanT<-Tempmean
meanP<-c()
for (i in 1:n) meanP[i]<-mean(c(Pressuremax[i],Pressuremin[i]), na.rm = T)
meanRH<-RHmean
#
#-------------------------------------------------------------------------------
# mean values
#-------------------------------------------------------------------------------
#
d<-3 # mean on d days
for (i in 1:(n-d)) {
  meanT[i]<-mean(meanT[i:(i+d)], na.rm = T)
  meanP[i]<-mean(meanP[i:(i+d)], na.rm = T)
  meanRH[i]<-mean(meanRH[i:(i+d)], na.rm = T)
  Y[i]<-mean(Y[i:(i+d)], na.rm = T)
}
AT<-273.15+meanT
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*T^4) 
#-------------------------------------------------------------------------------
#
T4_re<-lm(Y~I(AT^4))
summary(T4_re)
#
#-------------------------------------------------------------------------------
# Regression (model x1 + x2*t + x3T^4) 
#-------------------------------------------------------------------------------
#
t_T4_re<-lm(Y~meanT+I(AT^4))
summary(t_T4_re)
Advertisement
Una vita ___ interrotta.

Leave a Reply

Fill in your details below or click an icon to log in:

Gravatar
WordPress.com Logo

You are commenting using your WordPress.com account. ( Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. ( Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. ( Log Out /  Change )

Cancel

Connecting to %s