A simple example of mixed effects model using simulated data in R

First let's make the dataframe.

rm(list=ls())

library(lme4)

set.seed(1)

ID = c(rep('A',10), rep('B', 10), rep('c', 10))

x = c(1:10)

ya = x + 1 + rnorm(10)

yb = x + 2 + rnorm(10)

yc = x + 3 + rnorm(10)

df <- data.frame="" id="" x="rep(x,3)," y="c(ya,yb,yc))</font">

Then let's do liner regression.

mdl.lm <- span=""> lm(y~x, data=df)
coef(mdl.lm)

## (Intercept)           x
##    2.073487    1.001631

The linear regression cannot treat each ID seperately. It provides an overall intercept (2.073487) and overall slope (1.001631) for all the 30 data points.

Next, let's build a random-intercept mixed effects model.

mdl.mix.RI <- span=""> lmer(y~x+(1+1|ID),data=df, REML=F)
coef(mdl.mix.RI)

## $ID
##   (Intercept)        x
## A    1.287210 1.001631
## B    2.211162 1.001631
## c    2.722090 1.001631
##
## attr(,"class")
## [1] "coef.mer"

The above model can treat each ID differently, so the Intercept for each ID are different.

Last, let's build a random-intercept and random-slope mixed effects model.

mdl.mix.RIS <- span=""> lmer(y~x+(1+x|ID),data=df, REML=F)
coef(mdl.mix.RIS)

## $ID
##   (Intercept)         x
## A   0.9373545 1.0650834
## B   2.2293276 0.9929275
## c   3.0537802 0.9468823
##
## attr(,"class")
## [1] "coef.mer"

This model also treat each ID differently, and it gives different intercept and slope for each ID.

A GREAT way to reduce the .TIF file size while keep same high quality

Use Microsoft Office Picture Manager to open the original .TIF file, then export, set export size to 100%, DONE!

I had a >70MB .TIF image with 10000*10000 pixel, after applying this method, the image was shrinked to ~3MB. It's still 10000*10000 pixel and the quality is still good. That's amazing!

Use Sigmaplot to export high quality .TIF/.TIFF file

Sigmaplot is a powerful tool for drawing scientific charts, and it is also powerful at exporting high quality .TIF/.TIFF file, which can be directly submitted to the journal editors. Here's what I found out how to get high quality images from Sigmaplot:

Firstly, scale down the original graph. The default size of charts are relatively big, usually ~100mm*100mm. It's better to scale them down to 30%. Then select the items, and use highest DPI possible (for me, 600DPI). And change the exporting size to about 3 -5 time larger than the chart size. This will give you a large and clear TIF image.

AutoCAD TIF output quality Experiments

Since I don't really understand the relationship between DPI, PPI, image size, color mode, and the image quality, I decided to test them out by myself. Although it took me a whole nigh last night, I thought it is worth doing. Now I do know how to get a .TIF image with good quality (no blurry edges of letters and lines) and appropriate image size (several thousand pixels in each dimension, and the file size is about several tens MB).

Here's some tricks what I found out:

1. Scale down the drawing in AutoCAD. This won't decrease the image quality. Actually this can help decrease the file size of the output .TIF file. I used same settings for two prints, one is 10 times scaled up of the other one. And the two files are ~2MB and ~70MB respectively.

2. Use smaller print size in AutoCAD. At first I though I should print a very large size, maybe A2 or A1, to get a good image. Actually this has nothing to do with the image quality. When this size is large, the final .TIF image will be very large. I got a 10,000*10,000 pixel .TIF and my computer got stuck!

3. Usually 600dpi or 720dpi will be good enough.

4. Use 24Bit color. Other color mode doesn't work as well as 24Bit color.

A computer experiment on statistics: confidence interval for mean and variances of normal population

I am reading Statistics for Science and Engineering by Kinney today.In this book, there is a computer experiment assignment like this:

Select 1000 samples, each of size 10, from a N(10,5) distribution. Calculated the mean for each and a 95% confidence interval for u for each sample. Count the number of these confidence intervals that actually contain the true mean, 10.

I used MATLAB to do this.

clc;
clear all;
close all;
% Generate a normally distributed population.
Po=normrnd(10,5,[1,100000]);
SampleSize=10;
SampleNumber=10000;
%plot Po and histogram Po
figure
plot(Po);
figure
hist(Po, 100);

%take samples, each of size 'SampleSize'
Sa=[];
for i=1:SampleNumber
    for j=1:SampleSize
        index=round(abs(randn(1)*(length(Po)/10-1)))+1;
        Sa(i,j)=Po(index);
    end
end

%calculate the 95% confidence interval of each sample
mean=[];
interval=[];
for i=1:SampleNumber
    mean(i)=sum(Sa(i,:))/SampleSize;
    interval(i,1)=mean(i)-(1.96*5/sqrt(SampleSize));
    interval(i,2)=mean(i)+(1.96*5/sqrt(SampleSize));
end

%count the samples that contain 10
count=0;
for i=1:SampleNumber
    if interval(i,1)>=10 || interval(i,2)<=10
        count=count+1;
    end
end
count/SampleNumber

The results are around 0.05, which means that those confidence intervals have 5% chance not containing the true mean, 10. That's why those are 95% confidence intervals.


Plot of the population:

Histogram of the population:

Display message in Command Window, not using a pop-up message box

I like to use the input function very much. And I also like to give some instruction to other people who use my code, on what to input, in the command window, but not in a pop-up message box. After a little googling, I found the disp function.

clc;
clear all;
tic;
disp ('Hello, World!');
h=waitbar(0,'Please wait..');
n=0;
for i=1:100
    waitbar(i/100)
    for j=1:100
        for k=0:100;
            n=factorial(2);
        end
    end
end
close(h)
toc

MATLAB progress bar, show the progress of computing

Sometimes during a lengthy procedure, we don't have good way to determine if the code is still running or the computer got stuck. Using a progress bar will let you know approximately how long you have to wait til the run is over. Very cool!

clc;
clear all;
tic;
disp ('Hello, World!');
h=waitbar(0,'Please wait..');
n=0;
for i=1:100
    waitbar(i/100)
    for j=1:100
        for k=0:100;
            n=factorial(2);
        end
    end
end
close(h)
toc

xlswrite: output data to excel file {MATLAB functions}

Sometime when a column or a row of data is too long, you can't copy it from MATLAB and then paste directly into an Excel file. A much simpler way to do this is to use the xlswrite function. Just add a line at the end of your code:

xlswrite ('FileName.xls', VariableName)

If you want to put the data into sheet2 in the file, one more input argument is needed:

xlswrite ('FileName.xls', VariableName, 2)

Pretty easy, right?

Linear least-square regression

Of course, the linear least-square regression is an easy thing to do in Excel. Just adding a linear trend-line to the plot will do that for you. Anyway, I wrote this piece of code to do it in Matlab. 

% Single least-square linear regression
%Created on 02/17
x=input ('x='); %use [] to input a row
y=input ('y='); %use [] to input a row
X=mean(x);
Y=mean(y);
beta=(sum((x-X).*(y-Y)))/sum((x-X).^2);
alpa=Y-beta*X;
Equation=strcat('y=',num2str(beta), 'x+',num2str(alpa))

The num2str function was used to convert the beta and alpa, which are numbers, into strings.
The strcat function was used to concatenate the strings.

Calculate Sequence and Sum

For a(1)=sqrt(2) and a(n+1)=sqrt(2+sqrt(n)), calculate the sequence and sum for n=10.

clc; clear all; close all;

format long;

a=[];

S=[];

a(1)=sqrt(2);

for k=2:11

a(k)=sqrt(2+a(k-1))

end

a=a';

for m=1:k

S(m)=sum(a(1:m));

end

xlsread

Today I wrote a small piece to read data from xls files. Each file has several sheets and I only need one single column from each sheet. And at the end, I concatenated the columns to make a very very long colmn.

clc; clear all; close all;
O51=xlsread('TX2005.xls',1,'E1:E8760');
O52=xlsread('TX2005.xls',2,'E1:E8760');
O53=xlsread('TX2005.xls',3,'E1:E8760');
O54=xlsread('TX2005.xls',4,'E1:E8760');
O55=xlsread('TX2005.xls',5,'E1:E8760');
% O56=xlsread('NJ2005.xls',6,'E1:E8760');
% O57=xlsread('NJ2005.xls',7,'E1:E8760');
% O58=xlsread('NJ2005.xls',8,'E1:E8760');
O61=xlsread('TX2006.xls',1,'E1:E8760');
O62=xlsread('TX2006.xls',2,'E1:E8760');
O63=xlsread('TX2006.xls',3,'E1:E8760');
O64=xlsread('TX2006.xls',4,'E1:E8760');
O65=xlsread('TX2006.xls',5,'E1:E8760');
O66=xlsread('TX2006.xls',6,'E1:E8760');
O67=xlsread('TX2006.xls',7,'E1:E8760');
O68=xlsread('TX2006.xls',7,'E1:E8760');
O69=xlsread('TX2006.xls',7,'E1:E8760');
O71=xlsread('TX2007.xls',1,'E1:E8760');
O72=xlsread('TX2007.xls',2,'E1:E8760');
O73=xlsread('TX2007.xls',3,'E1:E8760');
O74=xlsread('TX2007.xls',4,'E1:E8760');
O75=xlsread('TX2007.xls',5,'E1:E8760');
O76=xlsread('TX2007.xls',6,'E1:E8760');
O77=xlsread('TX2007.xls',7,'E1:E8760');
O78=xlsread('TX2007.xls',4,'E1:E8760');
O79=xlsread('TX2007.xls',5,'E1:E8760');
O710=xlsread('TX2007.xls',6,'E1:E8760');
O711=xlsread('TX2007.xls',7,'E1:E8760');
O712=xlsread('TX2007.xls',7,'E1:E8760');
TXHourly=[O51;O52;O53;O54;O55;O61;O62;O63;O64;O65;O66;O67;O68;O69;O71;O72;O73;O74;O75;O76;O77;O78;O79;O710;O711;O712];

2nd version of the epa ozone data code

 %This program runs faster than the 1st version
%The data output is in a better form, no cell arrays are involved.
%When the data record is not a number but a string, it converts that record into NaN. In the 1st version, it onverts that into 0.

%This program loads a series of file of ozone data
clear all; close all; clc;
[TXSt TXCo TXSi]=textread('TXSites.txt','%n %n %n','delimiter','|');
[NJSt NJCo NJSi]=textread('NJSites.txt','%n %n %n','delimiter','|');
[LASt LACo LASi]=textread('LASites.txt','%n %n %n','delimiter','|');
FN=textread('filenumber.txt','%s');
TXData=[];NJData=[];LAData=[];%define arrays
for i=1: 161
    i
    [St Co Si Da Le]=textread(FN{i},'%*s %*s %n %n %n %*s %*s %*s %*s %*s %n %*s %s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s','delimiter','|');
    [r c]=size(St);
for j=1:r
   if (ismember(St(j), TXSt))&&(ismember(Co(j), TXCo))&&(ismember(Si(j),TXSi))
       TXData=[TXData; [St(j), Co(j), Si(j), Da(j), str2double(Le(j))*1000]];
   else if (ismember(St(j), NJSt))&&(ismember(Co(j), NJCo))&&(ismember(Si(j),NJSi))
           NJData=[NJData; [St(j), Co(j), Si(j), Da(j), str2double(Le(j))*1000]];
       else if (ismember(St(j), LASt))&&(ismember(Co(j), LACo))&&(ismember(Si(j),LASi))
               LAData=[LAData; [St(j), Co(j), Si(j), Da(j), str2double(Le(j))*1000]];
           end
       end
   end
end
end

The first code used to fetch the epa ozone data

%This program loads a series of file of ozone data
clear all; close all; clc;
TXSites=textread('TXSites.txt','%s');
NJSites=textread('NJSites.txt','%s');
LASites=textread('LASites.txt','%s');
year='2006';
NofF='161'; %number of total files
US='_'; %underscore
LQ='('; %lefrquote
RQ=')'; %rightquote
SF='.txt'; %surfix
n=textread('number.txt','%s'); %read 001, 002,003 ...
FN=strcat(year,US,n,LQ,NofF,RQ,SF);%generate FileName
State=[];County=[];Site=[];Date=[];Time=[];Level=[]; %define arrays
TXData=[];NJData=[];LAData=[];%define arrays
for i=1: 161;
FID=fopen(char(FN(i)));
A=textscan(FID,'%*s %*s %s %s %s %*s %*s %*s %*s %*s %n %s %n %*s %*s %*s %*s %*s %*s %*s %*s %*s %*s... %*s %*s %*s %*s','delimiter','|');
State=A{1,1};
County=A{1,2};
Site=A{1,3};
Date=A{1,4};
Time=A{1,5};
Level=A{1,6};
[r c]=size(State);
 for j=1:r
    ID=strcat(State(j),County(j),Site(j));
   if ismember(ID, TXSites)
       TXData=[TXData; [State(j), County(j), Site(j), Date(j), Time(j), Level(j)]];
   else if ismember(ID, NJSites)
           NJData=[NJData; [State(j), County(j), Site(j), Date(j), Time(j), Level(j)]];
       else if ismember(ID, LASites)
               LAData=[LAData; [State(j), County(j), Site(j), Date(j), Time(j), Level(j)]];
           end
       end
   end
 end
end
beep;

Use Matlab to edit sounds, images and videos

imread and imwrite

wavread and wavwrite

video? I don't know. I am a newbie!

Matlab Basics : dealing with arrays

Matlab is so powerful at handling arrays. It can generate lots of different kind of arrays, take one or some of specific elements in array and change the value of the elements. Since the data and files saved in computer are all basically 0 and 1, I guess that every file in the computer can be interpreted as certain kind of array, and thus can be accessed and edited by Matlab. (Of course, it's my GUESS.)

How to generate array using Matlab?

mauanlly input: A1=[1 NaN 3; 4 NaN 6; 7 NaN 9] % NaN means Not a Number.
use functions:

A2=magic(3) % magic array
A3=zeros(3) % all elements are0
A3=ones(3) % all elements are 1