title: Help me understand Bayesian updating

text: Suppose I have 5 possible events which either happened or did not happen in the preceding 10 time periods. How do I figure out how probable any event is in the 11th period?                   t11 t10 t9  t8  t7  t6  t5  t4  t3  t2  t1     A       1   1   1   1   1   1   1   1   1   1     B       0   0   1   1   0   1   1   1   1   1     C       0   0   0   0   0   0   0   0   1   1     D       1   0   1   1   0   1   0   1   0   1     E       0   1   1   1   1   1   0   0   0   0


title: Holt-Winters and importance of R-square

text: Is R-square an important measure in Holt-Winters method?


title: Reducing number of levels of unordered categorical predictor variable

text: I want to train a classifier, say SVM, or random forest, or any other classifier. One of the features in the dataset is a categorical variable with 1000 levels. What is the best way to reduce the number of levels in this variable. In R there is a function called `combine.levels()` in the **Hmisc** package, which combines infrequent levels, but I was looking for other suggestions.


title: Common statistical methods for analyzing weather station data

text: I have been tasked with analyzing temperature time series from multiple weather stations. I have minute and hour resolution temperature data spanning three years at 4 stations and 35 stations respectively. The stations reside on ~1000 km^2 (southern tip of Vancouver Island) so they cover a relatively concentrated area. I am looking for additional "statistical" techniques I can apply to the data. You might ask: what exactly do you want to glean from your data? Well, I am a complete amateur in weather and time series analysis. Pretty much anything that gives useful interpretation of data. So far I have performed the following:   * Power spectral density on minute data   * Probability distribution plots on minute data   * Auto lag correlation and cross-correlation between stations on minute data   * Principal component analysis on hourly data Techniques I have come across:   * Canonical correlation analysis   * Cluster analysis   * ? Any suggestions for books/papers/websites explaining analysis and practical applications are also much appreciated. I have found "Statistical Methods in the Atmospheric Sciences" (2006) by Wilks to be quite informative. Thanks!


title: Anomaly Detection in a set of points

text: I have a set of points in a matrix of size 100 x 100(total 10000 points). I know that there are roughly 500 anomaly points in it. There is a corresponding truth file which contains the true anomalous points which is not available while building the algorithm. The goal is to maximize the F-measure of the anomalies returned. How do I approach this problem? The dataset of all points


title: Stock & Watson Dynamic Ordinary Least Squares in R

text: I ran across the Stock & Watson 1993 paper on DOLS and was wondering if an implementation of it exists in R. A colleague was showing me EViews the other day and they have a DOLS model built in as part of their base package, but I have been unable to find a similar capability in any existing R package. Is there a package in R that implements the algorithm or is there a way to manipulate dynlm formulas to achieve the same effect?


title: calculate VIF for glmm

text: If we have colinearity among our axplanatory variables we can calculate a varianvce inflation factor to estimate the effects on our standard errors $$ VIF = 1/(1-R^2) $$ I am unsure how to calcuate this for a mixed effects model. I found this link but it is somewhat opaque (to me). Can someone provide me with a worded explanation of how to calcuate the VIF for mixed effects models?


title: Computing confidence region for Gaussian mixture model

text: I have a 2-d Gaussian mixture model and would like to compute a confidence region for it. Our application is that the two dimensions are latitude and longitude; that is, we want to say something like "the model is 95% confident that the true point lies within this (perhaps noncontiguous) region". In other words, I think what we want is to find the contour which contains 95% (or whatever) of the probability volume in the model. My questions:   1. Is there a simple and fast way to do this?   2. Are there any libraries which can do this, ideally in Python? We need the region as a (multi)polygon object which we can then pass around the rest of our software, not a plot.


title: Sensible to include ratio as a variable in logistic regression?

text: I'm creating a generalised linear regression using a binomial link function for two variables A and B. From looking at the data it appears that A/B may have discriminatory effect. Is it sensible to include A/B as an additional term in the model? Thinking about this myself, this appear term would be the inverse of including the interaction term A:B and therefore not the same as including the interaction. **EDIT 1:** Having thought about this a little bit more, if I include variables on a log scale would this be the same as including a ratio when the parameter estimates are negative? **EDIT 2:** I am experimenting with machine learning and have been exploring example data sets in the R package. I am currently using the `data(crabs)` in the MASS package and trying to predict crab gender and species from morphological variables. I first looked at the PCA of the variables to identify underlying vectors within the data. The first component in the PCA is all positive with respect to each variable, and from my understanding of PCA this means the variables are all correlated along this component. The second component contains both positive and negative values, specifically carapace width and rear width. A dot plot and density plot of these variable indicates two possible distributions. ![Plot of second component](http://bioinformatics- zen.s3.amazonaws.com/principal_components_analysis/second_component_dotplot.png) ![Plot of second component density](http://bioinformatics- zen.s3.amazonaws.com/principal_components_analysis/second_component_density.png) Is it therefore sensible to use the following logistic model for predicting species:               species ~ CW + RW + I(CW/RW)      As Stuart mentioned in the comments, these variables may be correlated with each other and therefore non-independent. I used bootstrapping of the crab data to test the models' accuracy. For the case of predicting sex, including the ratio term does seem to improve the model's predictive accuracy slightly. ![Model Bootstrapping](http://bioinformatics- zen.s3.amazonaws.com/principal_components_analysis/accuracy.png) I should also note that including all available measured variables improves the model accuracy but I would like to try and "infer" what the important variables may be using PCA.


title: Issue with Beta distribution

text: > **_Question 1_** :   >  When we merge two beta distribution, how do we arrive at the posterior > parameter values? Which is the right formula and why?   > In one of my earlier queries, Henry mentioned **Option 2** as the correct > way of calculating the posterior parameters, but I have seen a number of > papers giving **Option 1**. Hence, very confused. > >> **Option 1** : $\alpha (posterior) = \alpha (prior) + \alpha (likelihood) -1 $ or,   >  **Option 2** : $\alpha (posterior) = \alpha (prior) + \alpha (likelihood) > $. > > **_Question 2_** :   >  (Assuming the correct answer of **_Question 1_** is **Option 1**.   >  If we have the the prior distribution as $\beta (0.33, 1.35)$ and the > likelihood as $\beta (0.58, 0.88)$, then the posterior distribution is > $\beta (-0.09, 1.23)$. Hence, the posterior distribution does not follow a > beta distribution.   >  Is it the right way of combining beta distributions?   >  What do we do in this case? Please help.


title: Generalized estimating equations, correct use and interpretation (SPSS)

text: I am trying to predict a choice (which software will the participants use?) using the psychological measurements taken when the participants tried out the four possibilities. I have used linear mixed models, but the DVs are not suitable for that. GEE was suggested, but as I've never used it, I'm not sure it's proper, nor am I sure how to interpret the results. I read the most voted answers about GEE, but couldn't find answers. Details: The participants tried out four programs (in randomized order), and different responses were measured. So we have repeated measures and hierarchy with the four programs (lvl1) under the participants (lvl2), four rows per participant in the data. Right after the experiment the participants were given the initial choice: which ones would they use? Because the participants chose only one or two out of four possibilities, DV1 is a binary variable, this program chosen by this participant or not, for each row. Then, the participants were given the four programs and their actual use was recorded (in seconds) for an extended period of time. DV2 is therefore a continuous variable with a very bad skew, a lot of zeroes and low numbers to only a few high numbers. What I want to know is how the different response measures predict the actual use. Normally I would do a linear mixed model with something like this (using SPSS):               MIXED DV2 BY program WITH order predictor1 predictor2     /FIXED= program order predictor1 predictor2 program*predictor1 program*predictor2 | SSTYPE(3)     /RANDOM= intercept | SUBJECT(participant) COVTYPE(VC)     /REPEATED=order | SUBJECT(participant) COVTYPE(CS).      Where, to my understanding, the fixed effect for program tells me whether there are differences in program usages per se, fixed effects for predictors whether they predict regardless of program, and fixed effects for interactions whether the predictors predict differently depending on the program. But DV1 is multinomial and DV2 is horribly skewed. I've seen this done by transforming or standardizing DV2, but that does not really fix the distribution (results in one peak at zero and other one higher). I was told that tweedie distribution is what I'm looking for with DV2, and hence, GEE, which can also handle multinomial distributions. So now I got this:               GENLIN DV2 BY program WITH order predictor1 predictor2      /MODEL program order predictor1 predictor2 program*predictor1 program*predictor2  intercept=YES DISTRIBUTION=TWEEDIE(1.5) LINK=LOG      /REPEATED SUBJECT=participant WITHINSUBJECT=order SORT=YES CORRTYPE=INDEPENDENT ADJUSTCORR=yes COVB=robust UPDATECORR=1      /MISSING CLASSMISSING=EXCLUDE      /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION (exponentiated) WORKINGCORR.      (and DISTRIBUTION=MULTINOMIAL for DV1) 1) Is this a proper analysis for the described problem? I don't have the random designation in GEE, but - apart from the distribution, which is the point - this should be the same. However, the answer here says that the level of inference is different, but I really don't understand that. Does this mean that I can't get the program-specific estimates (equivalent of the interaction in LMM) from GEE? What are the Bs for interactions shown in the parameter estimates table then? 2a) If GEE is a correct analysis here, how should I interpret the betas? Even the highest numbers are very low: with example predictors, the largest B is 0.159, exponentiated to laughably small 1.172, while the intercept is 9.6, exponentiated to the right ballpark with 14 848 (seconds during three weeks). There is a "(Scale)" in parenthesis at the bottom of this table, never explained in the SPSS help files, that is 120 in this case (DV2 - with DV1 the scale is 1 as in SPSS documentation examples). Is this some kind of multiplication factor or something that I should apply to the exponentiated Bs to get the estimates in the same units as the original variables? 2b) Also, when the B is negative, the exponentiated B obviously only goes smaller than 1 but never below 0. Does this mean that GEE (with log link?) can never predict a **negative** association between the predictor and the dependent? If so, this is a huge problem, because I don't have any reason to assume that a rise in a predictor could never predict smaller usage. 2c) With DV1, the Bs are log odds, but the numbers are very similar to Bs from analysis with DV2, even if the latter ones should not be log odds. Am I doing something wrong?


title: What are the differences between the linear regression and multiple linear regression?

text: I'm interested to know that, what is the difference between linear regression and multiple linear regression? both of them seems same to me.


title: How to apply a logarithmic trendline if X Values = 0?

text: OK, I know its not possible because Log(0) is undefined. But here is my Problem: I measured a Weight gain of my Test object's over Time. Test objects, made of nonwoven material, were placed in very humid Atmosphere (enclosed above a Tank of heated Water). As they absorbed kondensated Water their total Weight would increase. Later the increase in Weight was lesser than in the beginning since test objects became saturated. ![enter image description here](http://i.stack.imgur.com/WoyKa.png) I plotted a graph in excel with time as x axis and weight as y axis. All of my test objects followed a similiar pattern. In order to take measurements i had to remove them from this humid enviroment for a short time. But since water vaporized really quickly the outcome is distorted. One would expect a logarithmic weight gain, but in a short time objects started to absorb water so slowly that the weight gain was slower than the weight loss(vaporization). Therefore 3. and 4. Point should be slightly higher than they are now, had I not disturbed the Process with my measurements. This would result in a logarithmic weight gain. ![enter image description here](http://i.stack.imgur.com/w9mNm.png) I woudl like to explain it in my report and provide a Trendline, but since my first Measurment was made at `t = 0 h` excel cannot calculate it. Here is my question: What is the most scientifically correct way to generate such a Trendline? I see a few possible options:   1. Ignoring the first value (but i don't like that),   2. Adding 1 hour to whole experiment (Whatever happend at `t = 0h` would be presented as it happened at `t = 1 h` and so forth). Or maybe there is another way? Any advice deeply appreciated. P.S. I am writing my Thesis for Dipl. Ingenieur in Mechanical engineering in Germany. EDIT: I added some additional info about my experiment. EDIT2: I added pictures and added even some more info. I hope it's now clear enough;)


title: How do you compare methods when they produce different binary answers?

text: I am trying to compare two methods for evaluating diagnostis of menisci in the knee. One is ultrasonographic where we can have 6 different diagnostic outcomes. The old method that serves as the gold standard also yields 6 different surgical outcomes. I would like to see if there is a correlation between the diagnostic and surgical findings. For instance - do you always get diagnostic finding type 4 when the surgical diagnosis is type B? Here is a sample of my cross-classification table: ![](http://i.stack.imgur.com/9rx2N.jpg) I have looked at all the different slides we have had in school but I can't find how to compare binary dataset with multiple measurements. I have included two different kind of schedules that could possibly enlighten the case a little bit. ![enter image description here](http://i.stack.imgur.com/A0QM4.jpg)


title: Earth Mover's Distance (EMD) between two Gaussians

text: Is there a close formula (or some kind of bound) on the EMD between $x_1\sim N(\mu_1, \Sigma_1)$ and $x_2 \sim N(\mu_2, \Sigma_2)$?


title: Comparing Twitter Trends on a Monthly Basis?

text: I'm trying to compare the frequency of a keyword is mentioned on Twitter in a monthly basis (e.g.: how many times the keyword "MissWorld" appears on September, how many it appears on August, how many it appears on July, etc). Something similar with Google Trends, or TwitterCounter, but for a keyword (not username). Would be great if it's possible to choose the location of the tweets. Is there any free tool available to do this? Or at least, an API?


title: Extracting values from R output of bugs command

text: I asked this question on http://stackoverflow.com/, but I couldn't get what I want. So, I am asking it here. I am running winbugs from R and I need to use some variables in R output. When I type `schools.sim\$mean\$theta[1]` in R, I get 10.2. However, when I type `schools.sim\$2.5%\$theta[1]`an error message come up. Any one what I am doing wrong or any other way to get the bayesian intervals? here is an example: This is the R code               library(R2WinBUGS)      data(schools)     J <- nrow(schools)     y <- schools$estimate         sigma.y <- schools$sd     data <- list ("J", "y", "sigma.y")               inits <- function(){      list(theta = rnorm(J, 0, 100), mu.theta = rnorm(1, 0, 100),      sigma.theta = runif(1, 0, 100))      }           schools.sim <- bugs(data, inits, model.file = "D:/model.txt",      parameters = c("theta", "mu.theta", "sigma.theta"),      n.chains = 3, n.iter = 1000,      bugs.directory = "D:/PROGRAMLAR/WinBUGS14/")           schools.sim      and this is the winbugs code which must be stored as model.txt in D.                model {      for (j in 1:J)       {      y[j] ~ dnorm (theta[j], tau.y[j])      theta[j] ~ dnorm (mu.theta, tau.theta)      tau.y[j] <- pow(sigma.y[j], -2)      }      mu.theta ~ dnorm (0.0, 1.0E-6)      tau.theta <- pow(sigma.theta, -2)      sigma.theta ~ dunif (0, 1000)     }


title: How to test for serial correlation of a time series itself (not residuals)?

text: We know that Ljung-Box test can be used to test for the residuals of a fitted model. But to test for the serial correlation of a time series itself, is there a way to do that?


title: Individual differences scaling, additional investigations

text: I am learning individual differences scaling (AKA three-way multidimensional scaling) and I want to know what different ways there are to state that my results are reliable. I had to perform individual differences scaling in 2 dimensions on the data, and give the group object space. I work in `R`, with the `smacof` package (this is obligatory). So far I am ok. But this is an additional question: > It is recommended to supplement the analysis with additional investigations > before placing too much belief in its results. List 3 investigations that > you would undertake and explain the reasons for performing these analyses. The answer I can think of: > 1) Check the overall stress of the model against a similar model with random > data. If the stress levels are not significantly lower, you probably got a > configuration that is based on coincidence. > > 2) Try the model with more and less dimensionalities, and look at the scree > plot of the stress levels. This can tell you if choosing a 2-dimensional > solution was the best decision. > > 3) Check the individual configurations, to see whether some persons used > different criteria than others. Then it would be better to look at the > individual pictures, instead of the group object space. Could anyone give any feedback on this answer? I can't seem to find much info anywhere else.


title: Calculate a tennis players chance to win a 5set match knowing he has 50% chance of winning 1SET?

text: Is the following correct? My friend tells me i need to (add up the binomial probabilities of 2 wins in 2 sets, 2 wins in 3 sets and 2 wins in 4 sets THEN multiply the sum obtained by p) P(3) = 5C3 x 0.5^3 x 0.5^2


title: Transforming two normal random variables

text: I'm reviewing for a test, and I am not sure if I am getting the right solution. Let $X$ and $Y$ be iid $\mathcal{N}(0, \sigma^2)$ random variables. a. Find the distribution of $U = X^2 + Y^2$, $V = \frac{X}{\sqrt{X^2 + Y^2}}$, b. are $U,V$ independent? c. Suppose $\sin(\theta) = V$. Find distribution of $\theta$ when $0 \le \theta \le \pi/2$. (tentative answers): I get   1. $f_{U,V}(u,v) = \frac{1}{4\pi} \sigma^{-2} \exp \left[ -u/(2\sigma^2) \right]| (1-v^2)^{1/2} + (1-v^2)v^2|$,    2. yes (density factors and supports dont rely on each other)and    3. $g(\theta) = \left[\cos^2(\theta) + \cos^3(\theta)\sin^2(\theta)\right]\frac{1}{8 \pi \sigma^4}$.  Anybody recognize any of these distributions?


title: Volatility of x and y variables in linear regression?

text: I have a simple regression with price returns: $r_{t+1} = \alpha + \beta r_t + \epsilon$ My question is: do I need to do anything if $r_{t+1}$ and $r_t$ are over different horizons? Suppose the _x-variable_ is a one month return but the _y- variable_ is a 1-day return: do I need to adjust the variables to account for the different volatilities? _Edit:_ In response to the current set of questions: Imagine that we have two time series $r_{monthly}$ with rolling monthly returns (sampled daily) and $r_{daily}$ with daily returns. And then we have the simple regression: $r_{daily, t+1} = \alpha + \beta r_{monthly, t} + \epsilon$ So we are predicting the next day return with the preceding monthly return. These returns are of different volatilities (assuming returns are gaussian, by the square-root of time rule). Should I do a transformation to normalize the volatilities before running the regression?


title: How to create common effect size estimate for meta-analysis based on a mix of mean changes and raw group means and standard deviations?

text: I am trying to do a meta analysis, involving continuous data. I am trying to do this meta analysis in Stata using the `metan` command. However, this data is organized in different fashions, where I have some studies with data listed as change in means with their standard errors and other data which is listed as the raw means and standard deviations. Is there a way to combine this data and proceed with the analyses?


title: Convert SAS NLMIXED code for zero-inflated gamma regression to R

text: I'm trying to run a zero-inflated regression for a continuous response variable in R. I'm aware of a gamlss implementation, but I'd really like to try out this algorithm by Dale McLerran that is conceptually a bit more straightforward. Unfortunately, the code is in SAS and I'm not sure how to re- write it for something like nlme. Would very much appreciate your help! Mike The code is as follows:               proc nlmixed data=mydata;       parms b0_f=0 b1_f=0              b0_h=0 b1_h=0              log_theta=0;                 eta_f = b0_f + b1_f*x1 ;       p_yEQ0 = 1 / (1 + exp(-eta_f));                 eta_h = b0_h + b1_h*x1;       mu    = exp(eta_h);       theta = exp(log_theta);       r = mu/theta;                 if y=0 then          ll = log(p_yEQ0);       else          ll = log(1 - p_yEQ0)               - lgamma(theta) + (theta-1)*log(y) - theta*log(r) - y/r;                 model y ~ general(ll);       predict (1 - p_yEQ0)*mu out=expect_zig;       predict r out=shape;       estimate "scale" theta;     run;      From: http://listserv.uga.edu/cgi-bin/wa?A2=ind0805A&L=sas-l&P=R20779 **ADD:** Note: There are no mixed effects present here - only fixed. The advantage to this fitting is that (even though the coefficients are the same as if you separately fit a logistic regression to P(y=0) and a gamma error regression with log link to E(y | y>0)) you can estimate the combined function E(y) which includes the zeroes. One can predict this value in SAS (with a CI) using the line `predict (1 - p_yEQ0)*mu` . Further, one is able to write custom contrast statements to test the significance of predictor variables on E(y). For example, here is another version of the SAS code I have used:               proc nlmixed data=TestZIG;           parms b0_f=0 b1_f=0 b2_f=0 b3_f=0                 b0_h=0 b1_h=0 b2_h=0 b3_h=0                 log_theta=0;                       if gifts = 1 then x1=1; else x1 =0;             if gifts = 2 then x2=1; else x2 =0;             if gifts = 3 then x3=1; else x3 =0;                     eta_f = b0_f + b1_f*x1 + b2_f*x2 + b3_f*x3;           p_yEQ0 = 1 / (1 + exp(-eta_f));                eta_h = b0_h + b1_h*x1 + b2_h*x2 + b3_h*x3;           mu    = exp(eta_h);           theta = exp(log_theta);           r = mu/theta;                if amount=0 then              ll = log(p_yEQ0);           else              ll = log(1 - p_yEQ0)                   - lgamma(theta) + (theta-1)*log(amount) -                      theta*log(r) - amount/r;                model amount ~ general(ll);           predict (1 - p_yEQ0)*mu out=expect_zig;           estimate "scale" theta;         run;       Then to estimate "gift1" versus "gift2" (b1 versus b2) we can write this estimate statement:               estimate "gift1 versus gift 2"       (1-(1 / (1 + exp(-b0_f -b1_f))))*(exp(b0_h + b1_h)) - (1-(1 / (1 + exp(-b0_f -b2_f))))*(exp(b0_h + b2_h)) ;       Can R do this?


title: How to convert a frequency table into a vector of values?

text: Using R or Excel, what is the easiest way to convert a frequency table into a vector of values? E.g., How would you convert the following frequency table               Value   Frequency       1.      2       2.      1       3.      4       4.      2       5.      1      into the following vector?               1, 1, 2, 3, 3, 3, 3, 4, 4, 5


title: Help with deploying a model

text: The two main packages I use at work are Palisade Risk and SPSS Clementine - they are both quite old versions and I've been supplementing the ability to analyze properly at work with more modern abilities available at home. E .g.at home I've been testing RapidMiner 5, which I really like in principle but have experienced a few issues with large data sets whereby the computational time progresses for hours then crashes for yet unknown reasons. What I like about Palisade & Clementine is they seem quite robust in ability to handle large amounts of input data and more importantly the ability to easily deploy models, i.e. once we have a model(s) that works we just feed it raw data and the outputs are predicted. However, Palisade & SPSS are both very limited in modern stats methods and machine learning algorithms that are hugely behind methods available in R/RapidMiner (Eg cross validation to name one). I would really like to convince my company to invest in R/Rapidminer but the ease of deployment of models is really the crux of the problem. May I ask how easy or difficult it is to achieve model deployment in the following and ask any advice from fellow professionals? i) RapidMiner ii) R iii) other recommendations?


title: What is wrong with lagged regressor in time series regression?

text: I am trying to explain a time series with the help of other related series. I really get nice fits using a standard LM approach with NeweyWest VC matrix. The fit even increases drastically when I replace some of the explanatories by higher lags (lag 10) of these explanatories (quarterly series) to the mix. Is it ok, to use such high lags? The length of the series is about 120. What's the risk of using high lags? At least it seems uncommon to me to use it.


title: Ordered sampling of uniform variables

text: I need to sample starting positions in a DNA sequence of given length L. These starting positions represent starting points of functional domains in DNA sequences. They have a length of d. The starting positions should be distributed uniformly and must have a minimal distance of d between them (functional domains should not overlap). Positions higher then L-d+1 are not possible because functional domains should be complete. The approach of sampling iterativly a starting position i and procecede with sampling in uniformly in the interval [i+d, L-d+1] does not produce a uniform distribution. Do you know of any solution to this?


title: What is the probability that a normal distribution with infinite variance has a value greater than its mean?

text: I got asked something similar to this in interview today. The interviewer wanted to know what is the probability that an at-the-money option will end up in-the-money when volatility tends to infinity. I said 0% because the normal distributions that underly the Black-Scholes model and the random walk hypothesis will have infinite variance. And so I figured the probability of all values will be zero. My interviewer said the right answer is 50% because the normal distribution will still be symmetric and almost uniform. So when you integrate from mean to +infinity you get 50%. I am still not convinced with his reasoning. Who is right?


title: How to set up a two-way mixed ANOVA with a covariate  in SPSS?

text: I wanted to carry out a two-way ANCOVA for my data. However, SPSS isn't liking that I have only one IV. I have   * One IV: Group (3 levels, i.e. 2 experimental groups and a control group),   * One DV: outcome measure (2 levels, i.e. time 1 and time 2),   * One covariate (I have checked and this demographic continuous variable moderately correlates with the DV). Please could someone help as to whether I am missing something in terms of putting this into SPSS?


title: Total waiting time at a server with exponential arrival times

text: **Question:** In a certain system, a customer must first be served by server 1 and then by server 2. The service time at server $i$ is exponential with rate $\mu_i$, $i = 1, 2$. An arrival finding server 1 busy waits in line for that server. Upon completion of service at server 1, a customer either enters service with server 2 if that server is free or else remains with server 1 (blocking any other customer from entering service) until server 2 is free. Customers depart the system after served by server 2. Suppose that when you arrive there is one customer in the system and that customer is being served by server 1. What is the expected total time you spend in the system? **Working:** Define: \begin{align*} X & = \text{Total time in system} \\ W_1 & = \text{Waiting time at server 1} \\ W_2 & = \text{Waiting time at server 2} \\ S_1 & = \text{Service time at server 1} \\ S_2 & = \text{Service time at server 2} \\ \end{align*} Then $X = W_1 + S_1 + W_2 + S_2$. We need to find $E[X] = E[W_1] + E[S_1] + E[W_2] + E[S_2]$. Since at the beginning there is only one customer in the system and that customer is being served by server $1$, then the waiting time at server $1$ is the same as the service time for the customer that is currently being served by server $1$. Therefore, $W_1$ is exponentially distributed with parameter $\mu_1$ and hence $E[W_1] = \frac{1}{\mu_1}$. Since $S_1$ is exponentially distributed with $\mu_1$, then $E[S_1] = \frac{1}{\mu_1}$. The Law of Total Expectation can be used to find $E[W_2]$. Define $A_1$ as the event where server $1$ completes before server $2$ and $A_2$ as the event where server $1$ completes after server $2$ completes. Then, \begin{align*} E[W_2] = E[W_2|A_1]P(A_1) + E[W_2|A_2]P(A_2) \end{align*} Define an ``arrival'' to mean the completion of service (at either servers), then consider the merging of the Poisson process represented by the service time at server $1$ and server $2$. The merged process is also Poisson with rate $\mu_1 + \mu_2$. Hence, \begin{align*} P(A_1) = \frac{\mu_1}{\mu_1+\mu_2} \\ P(A_2) = \frac{\mu_2}{\mu_1+\mu_2} \end{align*} By the memoryless property of the exponential distribution, the expected waiting time at server 2 given that server 1 has completed its service before server 2 is still exponentially distributed with parameter $\mu_2$. Hence, $E[W_2|A_1] = \frac{1}{\mu_2}$. If server 1 completes its service after server 2, then the waiting time at server 2 is zero, hence $E[W_2|A_2] = 0$. Thus, \begin{align*} E[W_2] = \frac{1}{\mu_2}\left(\frac{\mu_1}{\mu_1+\mu_2} \right) + 0 \left(\frac{\mu_2}{\mu_1+\mu_2} \right) = \frac{1}{\mu_2}\left(\frac{\mu_1}{\mu_1+\mu_2} \right) \end{align*} The service time at server 2 is exponentially distributed with parameter $\mu_2$, hence $E[S_2] = \frac{1}{\mu_2}$. In total, \begin{align*} E[X] & = \frac{1}{\mu_1} + \frac{1}{\mu_1} + \frac{1}{\mu_2}\left(\frac{\mu_1}{\mu_1+\mu_2} \right) + \frac{1}{\mu_2} \\ & = \frac{2}{\mu_1} + \frac{1}{\mu_2} \left(\frac{\mu_1}{\mu_1+\mu_2} +1 \right) \end{align*} **Query** I want to confirm whether the arguments I used are correct. More specifically, I am not too sure whether my argument to derive $P(A_1)$ and $P(A_2)$ is valid, can someone please confirm whether what I've done is right? Or provide any other kinds of solution would be cool too :)


title: ML vs MAP estimation, when to use which?

text: ML = Maximum Liklihood MAP = Maximum a-posteriori ML is intuitive/naive in that it starts only with _the probability of observation given the parameter (i.e. the likelihood function)_ and tries to find the parameter _best accords with the observation_. But it take into no consideration the prior knowledge. MAP seems more reasonable because it does take into consideration the prior knowledge through the Bayes rule. Here is a related question, but the answer is not thorough. http://dsp.stackexchange.com/questions/13174/differences-using-maximum- likelihood-or-maximum-a-posteriori-for-deconvolution-d So, I think MAP is much better. Is that right? And when should I use which?


title: Simulation of ARIMA (1,1,0) series

text: I have fitted the ARIMA models to the original time series, and the best model is ARIMA(1,1,0). Now I want to simulate the series from that model. I wrote the simple AR(1) model, but I couldn't understand how to adjust the difference within the model ARI(1,1,0). The following R code for AR(1) series is:               phi= -0.7048                                      z=rep(0,100)                                      e=rnorm(n=100,0,0.345)                            cons=2.1                                          z[1]=4.1     for (i in 2:100) z[i]=cons+phi*z[i-1]+e[i]        plot(ts(Y))                      How do i include the difference term ARI(1,1) in above code. Any one help me in this regard.


title: Replicating a plot from the ROCR Website

text: I'm trying to reproduce this plot from the ROCR website: ![enter image description here](http://i.stack.imgur.com/0ufkA.png) I can get something similar as follows:               library(ROCR)     data(ROCR.simple)     pred <- prediction(ROCR.simple$predictions, ROCR.simple$labels)     plot( performance(pred, "acc"))      ![enter image description here](http://i.stack.imgur.com/U0uoP.png) But I can figure out what `ROCR` metric to use to get postive vs negative densities, and I also can't figure out how to combine 2 `ROCR` curves into one plot. I figure since this figure shows up in several different places in the `ROCR` documentation, the code to replicate it must exist somewhere, correct?


title: Where can I use Chebychev's inequality?

text: What are some real-world applications where I can leverage Chebychev's inqueality? All the examples I find are either related to coin tosses or some school scores related problems. Are there any slightly-more complex use cases where this inequality "saves" the day?


title: Basic Simpson's paradox

text: I have a question about something that my statistics teacher said about the following problem: * * * There are two hospitals named Mercy and Hope in your town. You must choose one of these in which to undergo an operation. You decide to base your decision on the success of their surgical teams. Fortunately, under the new health plan, the hospitals give data on the success of their operations, broken down into five broad categories of operations. Suppose you get the following data for the two hospitals:               Mercy Hospital     Type         A    B      C    D      E    All     Operations  359  1836   299   2086  149  4729     Successful  292  1449   179   434   13   2366     Hope Hospital      Type          A   B  C   D   E   All     Operations   88 514 222 86  45   955     Successful   70 391 113 12  2    588      You notice that, in all types of operations, Mercy has a higher success rate than Hope, yet Hope has the highest overall success rate. Which hospital would you choose and why (choose two answers)? A) Mercy; since I would go in for a specific operation, I want the hospital that has the best success rate for that operation. B) Hope; since they do fewer operations in all categories, they are not "operation-happy" like Mercy. C) Hope; this is an example of Simpson's paradox and we should always chose the "obvious" conclusion. D) Mercy; looking at column E, Mercy clearly does more difficult surgeries and so is probably a better hospital. E) Hope; it has the better overall success rate. F) Mercy; this is an example of Simpson's paradox and we should always chose the opposite of the "obvious" conclusion. * * * **My question isn't even about the occurrence of Simpson's paradox in this situation. My question is simply about the fact that my professor insists that A) and D) are the right answers instead of A) and F). He says,** > "Because the success rate is so low for Type E surgeries,we can conclude > that they are difficult and not just uncommon. Hence, Mercy probably has > better equipment/doctors when compared to Hope." I don't understand how he could imply on a statistical basis that he can tell that Mercy does "more difficult surgeries". It is obvious that Mercy has better success rate at type E surgeries, but why does that mean they do "more difficult surgeries". I think I am being screwed over by the wording of this problem and the professor isn't budging. Can someone please explain why I am wrong or how I can explain this to the professor?


title: Interpreting logistic regression coefficients with a regularization term

text: I understand the coefficients of a logistic equation can be interpreted as odd ratio. If a regularization term is added to control for over-fitting, how does this change the interpretation of the coefficients?


What is the relation between estimator and estimate?


title: Perceptron overfitting?

text: I'm trying to judge the performance of my perceptron linear discriminant. In one instance I'm training on a sample size of 150 and on another I'm training on a sample size of 1500. I test both of these on the same testing set. Yet, the peceptron that is trained on the larger set performs more poorly. I'm not sure how to describe this...could this be an example of overfitting??


title: How to calculate barometric pressure trend?

text: Greetings, I need to determine the trend for barometric pressure based on measurements taken at 15-minute intervals over a 3 hour period. A rise or fall trend occurs if the the pressure changes more than 0.1 kPa, and rapid rise or fall if pressure changes more than 0.2 kPa over a three hour period. Measurements are in kPa, ranging roughly from 97.00 to 104.00. Haven't used my college math in over thirty years - just can't seem to remember this. Any help in this area would be greatly appreciated. Thanks!


title: Compute confidence intervals for multinomial proportion

text: I need to compute the confidence intervals (CI) for a multinomial proposition. I am not a statistician but I can understand basic stuff... Looking around I have found this resource that explain how to do it: Glaz, J. and C.P. Sison. 1999. "Simultaneous confidence intervals for multinomial proportions". Journal of Statistical Planning and Inference 82:251-262. I was able to understand more or less the approach but I cannot understand a detail. To compute the CI as reported in formula (8) in the paper I need to use their Theorem 2.1 which depends on a truncated Poisson distribution and in particular I need to compute the central and factorial moments of such distribution (mu_i, sigma_i, and mu_n_i). I do not understand how I can compute these moments....without them I cannot compute the intervals.... Thank you so much for your help! G.


title: Is "not-overfitting" a utopian scenario?

text: We say a model overfits when classification error increases on the test data. The reason behind this is that the training data is not a representative of the distribution from which data is sampled. Can we **ever** say that a model with a given set of parameters does not overfit? My question stems from the fact that the test data is also not a representative of the distribution, so having low test error also does not mean that the model will not overfit. Instead, we can only claim that **with high probability** this model will not overfit. Is there a flaw in this line of reasoning?


title: Statistical test for uniform distribution

text: I have a sample of 5 numbers from known interval [0, 10]. Is 5 numbers is enough to make some conclusions about whether these numbers are drawn from uniform distribution or not?


title: Suppose $X_1,...,X_n$ are iid from $N(\mu,\sigma^2)$, how can I find $P(T(X)>200)$?

text: Suppose $X_1,...,X_n$ are iid from $N(\mu,\sigma^2)$ with $\mu$ unknown and $\sigma = 2$ known. Letting $T(X) = \sum{(X_i-\mu)^2}$ from i = 1 to 40, how can I use the central limit theorem to approximate $P(T(X)>200)$? I found that $E(T(X))=40$ and $Var(T(X))=160$ and $T(X)$ ~ Chi square with n=40 df. The below is the correct answer, but I do not know how this makes sense since my understanding of the central limit theorem is that the average of iid random variables is normal, but how am I utilizing the central limit below? Thanks!!!! ![enter image description here](http://i.stack.imgur.com/6HNWe.png)


title: Triangular distribution

text: I am working on a dataset that ranges between 0 - 100 and typically is centered around 50 for the most part. It seems like the dataset shows a triangular distribution given its bound between 0 - 100. I wasn't sure and just wanted to double check if this would be a triangular distribution? Can we compute z-score for this dataset? If not how would you recommend analyzing this dataset on a cross-sectional as well as time series basis?


title: How to describe and present the issue of perfect separation?

text: Folks who work with logistic regression are familiar with the issue of perfect separation: if you have a variable specific values of which are associated with only one of the two outcomes (say a binary $x$ such that all observations with $x=1$ have outcome = 1), the likelihood blows up, and the maximum likelihood estimates run off to infinity. `glm` in R may or may not handle that terribly well, as the perfect prediction error message can show up for reasons other than perfect prediction/separation. `logit` in Stata identifies such variables and problematic values, and discards them from the analysis. My question is different from what to do if you have a perfect separation. That I can handle by recoding my variables (they are all categorical, so I can simply combine categories), or with Firth version of logistic regression if I want to be fancy. Instead, I wonder what the common ways are of _describing_ this. I have a data set with circa 100 patients with about 50% proportion "positive", and some categories of the demographic variables produce this perfect prediction. Let's just say that all 7 green-eyed people have a "positive" outcome. This may be a small sample peculiarity that would disappear if I had sample size of 1000 and 70 green-eyed people, but it may be clinically meaningful, as in that larger sample I could have 60 out of 70 green-eyed people that would have a "positive" outcome with high odds ratios. So it's nice to say that I used a Bayesian or some other shrinkage method, but in describing how I got there, I would need to admit that I had perfect prediction/separation, and had to find a more sophisticated technique to get any results at all. What would be a good language to use here?


title: Which regression and correlation to use?

text: I have a database of patients with their heart rate before and after surgery. I would like to do a few analyses:   1. Define a threshold value for an optimal heart rate (i.e., above which there is a significant risk of post operative complications).   2. Test this hypothesis: Does increasing heart rate increase chances of post operative complications?   3. If no 2 is true, then describe the increase in risk of post operative outcomes as heart rate increases. For example, I want to say for every 10 beat per minute increase in heart rate, the relative risk of post operative complications increases by 10%. What type of regression or tests would be suitable for the analyses above?


title: How to interpret and code semantic differential scale into SPSS?

text: Is the semantic differential rating scale below considered a scale of intensity? This is a sample of 7 items contained in a group PQ. Legend: 1 and 7: (extremely), 2 and 6: (moderately), 3 and 5:(fairly), 4: (Neutral)                              1  2  3  4  5   6  7      people centric | O | O | O | O | O | O | O | Technical. Isolates | O | O | O | O | O | O | O | Connects.   Pleasant | O | O | O | O | O | O | O | Unpleasant. Original | O | O | O | O | O | O | O | Conventional. Simple | O | O | O | O | O | O | O | Complex. Ugly | O | O | O | O | O | O | O | Pretty. 1) If a participant chooses 1 and 7, how are they coded to be input into SPSS? 1 is extremely people centric and 7 is extremely Technical Centric. If we take the average of the 7 items to obtain a mean value to input into SPSS, then 1 and 7 will have unequal weight age. I intend to perform a parametric analysis (say ANOVA)? Thank you


title: Correlation between vegetation and erosion

text: I am wanting to find the correlation between two variables: a categorical vegetation variable (with categories 1-4) and a continuous variable of erosion (with values anywhere between 0.2 and 9). I think I would be required to use a non parametric test to find the correlation, but bit of a newby to stats so could anyone suggest the best test to use? I also wanted to know how sample size can affect these tests, ie with a sample size of 100 or 1000?


title: Can I add means together for a set of variables related to a construct, so as several constructs can be compared?

text: I have 8 sets of variables, taken from another study, each set related to either an intrinsic or extrinsic motivational construct (e.g. "challenge" or "recognition"). I want to see which construct is strongest over a number of participants. Can I add the mean of each variable (over participants) within a set of variables (the construct) to ascertain this …or, is this statistical madness? How do I compare one set of variables with another to see which is strongest? (I.e. I want to compare constructs so as I can see which motivational set is strongest for a group of students).


title: What can I do with a modified Likert-scale questionnaire?

text: We are using a Likert scale questionnaire and need to be able to identify the rate of the students. What is the best way we identify this rate? What statistical methods are typically used? We based our motivation questionnaire on the MSLQ (pdf) but we just have a 5-point Likert response format instead of 7 points. We also remove some of the items in the questionnaire which is not related to our study. We want to know if the students were being motivated to the strategy that we used. And the answer will be based on the questionnaire that we have given them. How can we show it through statistic?


title: Self-study question on T-Distribution

text: I am currently doing a True/False exercise. One of the questions as follows: > A marine drill instructor recoded the time in which each of 20 recruits > completed an obstacle course both before and after training. To test whether > any improvement occurred, the instructor would use a t-distribution with 38 > degrees of freedom. While this statement is obviously **False** due to the fact that the degree- of-freedom cannot be 38 since there are only 20 recruits, I was wondering if the following statements will make the statement true:   * degree of freedom should be 19   * it is assumed that a population standard deviation is not available, hence the use of the t-distribution Appreciate any advice please.


title: How to interpret ADF-test optimal lag model reduction?

text: I was looking at howto figure out the optimal lag order for an ADF-test when I came across this PPT presentation. There are two examples of a ADF-test model reduction, one on p23 and one on p36. I have two questions:   1. How do I come up with the F numbers of the reduction? For example, on p23 it says "F(1,487) = 0.92201 [0.3374]"? I understand what an F number is, but what calculation lies behind the 0.92201 number (and 0.3374)?   2. How do I interpret the numbers? For example, why is model 4 the preferred one on p23? I'm not a complete noob, but I certainly wouldn't mind a nicely laid out pedagogical answer. Thanks!


title: Support Vector Machines....||w||, support vector coefficients, etc

text: I'm working with LibSVM and Matlab but am running into some issues. Specifically, I think I'm misunderstanding some theory. I've read here and there about the support vector coefficients, but I don't understand them. Are the support vector coefficients the ||w|| to that support vector?? Specifically, my long term goal is to find the distance from a support vector (or any training/testing vector really) to the hyper plane. Can anyone please shed light on all this! Thanks


title: Question on Portfolio Risk & Covariance

text: Good evening, I was reading a stats textbook where I came across the following statement: > In order to reduce the portfolio risk of a portfolio involving two risky > investments, we should choose their covariance to be negative, rather than > positive I don't really get what it means. From what I understand: > Cov = 0, the variables are not related > > Cov = 0 to 1, the if the value of one variable increases, so will the other > > Cov = -1 to 0, the if the value of one variable increases, the other will > decrease Am I right to then say that why having a negative covariance will reduce risk is because when one investment falls, the other will react to counter? So in such cases, the profits will be always minimal but safe. In contrast, if the investments are positively correlated, then the portfolio will be highly risky, but highly profitable if one makes a profit.


title: Prediction of continuous variable using "bnlearn" package in R

text: I use bnlearn package in R to learn the structure of my Bayesian Network and its parameters. What I want to do is to "predict" the value of a node given the value of other nodes as evidence (obviously, with the exception of the node whose values we are predicting). I have continuous variables.               library(bnlearn)                       # Load the package in R     data(gaussian.test)     training.set = gaussian.test[1:4000, ] # This is training set to learn the parameters     test.set = gaussian.test[4001:4010, ]  # This is test set to give as evidence     res = hc(training.set)                 # learn BN structure on training set data      fitted = bn.fit(res, training.set)     # learning of parameters     pred = predict(fitted$C, test.set)     # predicts the value of node C given test set     table(pred, test.set[, "C"])           # compares the predicted value as original      Now, this code works well and gives a table in which you can see that predicted values for node C are **exactly** same as the original value of node C in test set. I do not understand the reason for that, could someone please explain it? I know, I am providing entire df of test set that already has the value of node C in it. But if I give the data of other columns, it gives an error. So, I tried an alternative of putting other values to 0.               test.set$C = 0                     # To not give the original value of node C as evidence     pred = predict(fitted$C, test.set) # predicts the value of node C given test set     table(pred, test.set[, "C"])       # compares the predicted value as original      Is this approach wrong? (Using "NA" is not allowed.)


title: Self-study question on CLT

text: Good evening everyone, I am currently doing a self-study on CLT and was working on an exercise which asks if the following statement is true or false > CLT guarantees that the population mean is normally distributed whenever > sample size is sufficiently large. I can understand that this statement is **False** , but I am pretty weak with the concept, thus will like to take this opportunity to validate my thoughts. > CLT does not guarantee that the population mean is normal, instead it allows > us to treat a population as normal if a given population size is large > enough. Also, the size of the sample will not affect the actual distribution > of the population, but rather simply gives us a more accurate approximation. Appreciate if seniors and exports here can give me some guidance.


title: Methods for testing a Bayesian method's software implementation

text: I'm a frequentist but have been exploring Bayesian techniques lately. One thing that bothers me is that I'm never completely sure that I've implemented them correctly. To be specific, I'm not talking about convergence, burn-in time, or autocorrelation. I'm talking about, 'did I code the formulas correctly?'. What I do now is to jack up the sample size and watch for the posterior to hug the true parameter values (with simulated data). But that's completely ignoring the power of the prior and the great small sample properties that Bayesians love to talk about. Frequentist techniques make testable predictions in the form of relative frequencies. But in the Bayesian context, we lose the relative frequency-based interpretation of probability. So how does one programmatically test a Bayesian method's software implementation?


title: Sentimental Analysis using Naive Bayes

text: I am working on problem solution where I am collecting social feeds from twitter and Facebook for a product X . I am labeling these posts,comments or tweets using five labels               --Positive     --Negative     --Campaign      --Reply     --Queries       I have a training set of around 5000 which includes tweets,Facebook posts and Comments . But these training set are unbalanced and have more of Negatives and Campaign data . Below is the list of sentiments and their count:               --Positive--> 492     --Negative--> 2193     --Campaign--> 1422     --Reply--> 430     --Queries--> 922      I am using Naive Bayes for predicting these sentiments . As you can see the above training set is high unbalanced is there any way that I can improve my model with these training set . Any suggestion for improving my prediction model would be helpful . I am using Mahout for these building this prediction model . Thank you


title: How to test if the errors of a linear regression are constant?

text: > **Possible Duplicate:**   >  Kruskal-Wallis or Fligner test to check homogeneity of variances? Is there a test to check if the residuls of a linear regression are constant? Thanks


title: Multicollinearity of Socioeconomic index vs. wealth variables separately?

text: I am running Cox regressions in a large dataset to determine effects of individual and area characteristics on mortality outcomes. (My model is not a multilevel model, as the initial significant variation (random intercepts) between areas disappeared once SES variables have been adjusted for). I have three categorical individual level socio-economic measures (car ownership, housing tenure, and education). These can be used to construct a socioeconomic index with a scale of 0-7 (based on their categories). On the area level, I have a "morbidity" variable, which is the proportion of adults with an illness in the area. When using the socioeconomic variables on their own—i.e. not as an index—they are not correlated with each other, nor with area illness. However, the individual SES index aggregated on area level is highly correlated to area illness (0.78). Can I put in the three SES measures separately (not as an index) and the area level measure for illness in the same model, or is that problematic? Many thanks! Edit: I share your view that multilevel models are important to estimate the variance. However, a few things make me think I it might be ok to go with a non-MLM model in this case. To give a bit more info on my analysis: I have approximately 500 000 individuals clustered in 571 wards. I’m examining a variety of mortality outcomes, some of which are quite rare, which means I do not have a death/event in every ward during my follow up period. I think I need at least one event in every level 2 unit (in this case ward) unit to run a MLM (or Cox survival analysis with shared frailty)? Due to the large sample size Cox models with shared frailty (in Stata 12) take a very long time to run. Run instead with Streg my theta is 0.008 (st err 0.01) when all individual level factors (SES, age, family etc) are controlled for. I relatively new to survival analyses with shared frailty and haven’t come across a good explanation on how to interpret theta. IS it the same as level 2 variance estimated e.g. in a normal multi-level model in MlwiN or R? Any help here appreciated! I interpret this as when demographic and socioeconomic factors explain any variation between wards in risk of death. AIC for models without shared frailty (by ward) are lower than models without (of at least 2, often more). With this information in mind, would you still advice me to run models with shared frailty? Thank you again!


title: time series analysis: metrics which ignores length

text: i'm quite new here! I have a problem - set of time series with different lengths each. And i need to measure some similarity among them through metrics. I tried to use DTW metric, but it seems to take into account length. Please, could you give me some keys about metric, which i need to use? Thank you in advance! Any help would be nice.


title: Comparing 2 sets of observations

text: I am analysing signals from 5 sensors and comparing what they pick up between day and night. Each sensor records a time series and this is compared to a reference time series to determine whether what the sensor is picking up is significantly different to the reference signal. This is done a number of times (represented by Obs 1, Obs 2 etc.) This information is then represented in a matrix (as below) showing p values for each comparison. ## _DAY TIME_                         Obs 1    Obs 2    Obs 3    ...     Sensor A   0.01     0.11     0.01     Sensor B   0.10     0.23     0.91     Sensor C   0.4      0.31     0.56     Sensor D   0.005    0.001    0.002     Sensor E   0.02     0.002    0.011      ## NIGHT TIME                         Obs 1    Obs 2    Obs 3    ...     Sensor A   0.31     0.18     0.02     Sensor B   0.14     0.02     0.09     Sensor C   0.41     0.03     0.16     Sensor D   0.11     0.12     0.23     Sensor E   0.33     0.21     0.11      My aim is to compare Day Time and Night Time readings for each sensor and determine whether there is a difference between Day and Night. What is the best way to make this statistical comparison given that my observations are represented in the tables as p values?


title: Cluster-robust Standard Errors

text: acoording to paper http://cameron.econ.ucdavis.edu/research/Cameron_Miller_Cluster_Robust_October152013.pdf on page 4 the “failure to control for within-cluster error correlation can lead to very misleadingly small standard errors, and consequent misleadingly narrow confidence intervals, large t-statistics and low p-values”. However, unfortunately my cluster-robust standard errors not always get larger than the default standard errors; sometimes they become smaller. What's the reason for that?! Since cluster-robust standard errors are also heteroskedastic-robust, I was wondering if maybe these two issues outweigh each other? However, this would be weird as well since tests for heteroskedasticity were negative. So, does somebody has an explanation? Hope you guys have an answer. Regards, Phil


title: How to transform a data set to downweight (relative) small numbers

text: I've got a data set looking at how different groups change over time.               #High Abundance, Low Change     HALC<-c(100,99,101,98,99,100,100,101,99,100)          #Low Anundance, Low Change     LALC<-c(1,2,1,2,2,2,1,2,1,2)          #High Abundance, High Absolute Change     HAHAC<-c(100,99,98,91,86,50,45,30,21,9)          #High Abundance, High Absolute Change     LAHAC<-c(9,21,30,45,50,86,91,98,99,100)          #High Abundance, Changes 10%     HA10C<-c(100,101,102,103,104,104,106,107,108,110)          #Low Abundance, Changes 10%     LA10C<-c(1,1.01,1.02,1.03,1.04,1.04,1.06,1.07,1.08,1.1)          #High Abundance, Changes 100%     HA100C<-c(100,110,120,130,140,140,170,175,187,200)          #Low Abundance, Changes 100%     LA100C<-c(1.00,1.10,1.20,1.30,1.40,1.40,1.70,1.75,1.87,2.00)     DF<-c(HALC,LALC,HAHAC,LAHAC,HA10C,LA10C,HA100C,LA100C)               DF<-data.frame(HALC,LALC,HAHAC,LAHAC,HA10C,LA10C,HA100C,LA100C)     row.names(DF)<-c(1,2,3,4,5,6,7,8,9,10)      I'm running a secondary analysis that turns out to be much more sensitive to relative change within a single group, as opposed to the entire data set. For example, looking at these two groups:               #High Abundance, Changes 10%     HA10C<-c(100,101,102,103,104,104,106,107,108,110)          #Low Abundance, Changes 10%     LA10C<-c(1,1.01,1.02,1.03,1.04,1.04,1.06,1.07,1.08,1.1)      If I run my analysis on the above data, these two groups get weighted exactly the same because both change only 10%. But in terms of my interest in the dataset, I'm much more concerned with change in units (change relative to the entire data set). AKA That the High Abundance/ 10% change group increased by 10 units is much more important than the Low Abundance/ 10% change which only increased 1 unit. So can anyone recommend a data transformation that would either: 1 Transform the data to reflect the importance of change in relation to the relative size of the dataset (AKA Even though an increase from 1 to 2 is a 100% change, in relation to the data set it is a very small change) 2 or at least downweight the impact of low relatively small data Because in the end my secondary analysis is going to give more weight to something that increases from 1 to 2 than 100 to 102 (based on percent change), which I would like to avoid. As a final note, I need to avoid negative numbers and zeros in my final transformation. And still preserve the change over time aspect of the data (increasing or deceasing over time).


title: Good algorithms for feature extraction from images?

text: I am searching for some algorithms for feature extraction from images which I want to classify using machine learning . I have heard only about [scale- invariant feature transform][1] (SIFT), I have images of buildings and flowers to classify . Other than SIFT what are some good algorithms .


title: Help me prove that my pickup technique for women is statistically valid

text: My girlfriend and I were having a debate about a magic trick that could be used to pick up women on the street. I am a magician and I have a coin. I walk up to a complete stranger with the coin in my hand and know which side is up. I ask her to state heads or tails. If she picks the correct one I tell her, "I knew you would pick that one. You should marry me!". If however, she doesn't pick the one I know that I have in my hand facing up. I flip the coin, and if I get it correct, I instead say, "You should go and get it drink with me immediately!" What there chances that she will be told either that she should marry me or go out and get a drink with me? I have a 50% chance of her getting the first line 100% of the time. That leaves another 50% unaccounted for. I now have a 50% chance of her getting the second line, but only 25% of a chance over all. Can these two numbers be added together to get 75% chance total? My girlfriend says it doesn't work like that I am wrong.


title: Sampling from the conditional distribution assuming sampling from the joint

text: I am struggling with this question, which I thought it should be easy: suppose we have a method of sampling from the joint distribution of a collection of (discrete ordinal) random variables. We do not know the density. We need to find a way to sample from the **conditional** distribution of a subset of them, $X_i$, given specific values for the others, $Y_i = y_i$. Since they are discrete I thought of generating a very large sample from the joint and keep only those where $Y_i =y_i$. I am not sure if this makes sense and it also looks potentially extremely inefficient. Is there any other way to do this?


title: One component in PCA is always the mean vector in two-dimensions but not three

text: I've been testing PCA via SVD to decompose a simple time series data matrix, $X$. I have two signals $x_1(t)$ and $x_2(t)$ in a data matrix where $M$ rows represents each timepoint sample and each column represents $x_1$ and $x_2$. The mean signal, $\hat{x}$, is defined as the mean along the row axis (average of $x_1$ and $x_2$ along each timepoint). I normalize each column of $X$ by subtracting its mean and dividing by the standard deviation. When I use [U S V] = svd(Xz) in matlab, regardless of how the variables are distributed (whether they are correlated or uncorrelated), one of the columns of the right singular matrix, V, always points in the same direction (to a multiplicative constant) as the mean vector $(1/2, 1/2)$. But when add an additional third time series, this is never the case (where the mean vector is $(1/3, 1/3, 1/3)$. Because I normalize the standard deviation for each vector, it does make sense that the direction of most variance given by PCA would be the diagonal 45 degree line. But if both variables $x_1$ and $x_2$ are independent gaussians, couldn't the PCA direction be any direction since the distribution is radially symmetric? MATLAB Code:               s = RandStream('mcg16807', 'Seed', 0);     RandStream.setDefaultStream(s);          G = zeros(1000,2);     G(:,1) = 40*randn(1000,1)-100;     G(:,2) = tan(G(:,1)) + randn(1000,1);     X = G - repmat(mean(G),[size(G,1), 1]);     Xz = X./repmat(std(G),[size(G,1), 1]);     [U S V] = svd(Xz);     G_mean = mean(Xz,2);     corrcoef(G_mean, U(:,1))     corrcoef(G_mean, U(:,2))


title: Can the sample size of subgroups influence significance of correlations?

text: I am currently helping a friend with some statistical analysis. There is a main hypothesis, for which a significant relationship can be shown. Lets call the predictor variable A and the predicted variable B. Her sample size is 200 people. However the distribution of the sexes is very uneven, i.e. there are 150 women and only 50 men. There also is a significant relationship in her data between gender and the predicted variable B. Furthermore if all factors are analyzed together, the significance of the relationship between A and B drops below the alpha-level. Is it possible that the unevenness of the sample can somehow influence the result in this case? Or is there some other effect that explains this situation. If so what is this effect called (assuming it is known). Maybe the stronger relationship between gender and B is shadowing the relationship between A and B or something like that.


title: How to calculate the signal-to-noise ratio for an image?

text: I've started with a “perfect” image 'without noise'. Then I added some noise to it. I've used gaussian white noise for this purpose.               In = In+0.5*randn(size(Io))      The combined image (my “noisy image”) has a signal to noise ratio with some meaning because I can compare it to the perfect image. My question is how can I calculate the 'SNR Value'.


title: Multi-sample nonparametric test for matched data

text: I've been wracking my brain trying to come up with a test to do something similar to a paired wilcoxon rank sum test for more than two samples. I've explored the Friedman and Kruskal-Wallis tests, but neither of those do quite what I want (I think). The particular case in the data set I'm dealing with is the following. I have two (or more) samples with measurements that are matched (i.e., the first measurement from sample A corresponds to the first measurement from sample B). If I only had two samples, I could simply take the difference between each measurement (e.g., A(1) - B(1), A(2) - B(2) etc.) and use a Wilcoxon test to see if the observations are significantly different from 0. Does anyone have any suggestions on the most appropriate way to extend this kind of logic to more than 2 samples? Kruskal-Wallis does not take matching observations into account (as far as I can tell) and the Friedman test wouldn't detect a difference in the example below (which I would want to call different):               +-----------+-----+-----+-----+-----+     |Observation|1    |2    |3    |4    |     +-----------+-----+-----+-----+-----+     |Sample A   |10   |10   |0    |0    |     +-----------+-----+-----+-----+-----+     |Sample B   |0    |0    |10   |10   |     +-----------+-----+-----+-----+-----+      My thought is that given N samples I can take every possible pairwise difference between the samples, sum these differences, and use a wilcoxon test to compare the sum to 0. My naive statistical intuition is that since I'm taking differences and sums of values from the same distribution (the assumption in the wilcoxon test) these operations will result in a distribution that is the same as those of each of the samples (and consequently are valid). If any clarifications would help answer my question please let me know and feel free to point out if I'm completely off-base :). Thanks in advance for any advice.


title: Making big, smart(er) bets

text: I've been trying to code an algorithm to suggest bets in 1X2 (weighted) games. Basically, each game has a set of matches (home vs away teams):   * `1`: home wins   * `X`: draw   * `2`: away wins ![BWin 1X2 Betting Game](http://i.stack.imgur.com/mD1w8.jpg) For each match and symbol (`1`, `X` and `2`), I will assign a percentage that represents the chances / likelihood of that symbol being the correct match outcome. Here is an array representing the structure:               $game = array     (         'match #1' => array // stdev = 0.0471         (             '1' => 0.3,     // 30%     home wins             'X' => 0.4,     // 40%     draw             '2' => 0.3,     // 30%     away wins         ),              'match #2' => array // stdev = 0.4714         (             '1' => 0.0,     //   0%    home wins             'X' => 0.0,     //   0%    draw             '2' => 1.0,     // 100%    away wins         ),              'match #3' => array // stdev = 0.4027         (             '1' => 0.1,     //  10%    home wins             'X' => 0.0,     //   0%    draw             '2' => 0.9,     //  90%    away wins         ),     );      I also calculate the standard deviation for each bet (commented in the above snippet); higher standard deviations represent a higher certainty, while the matches with the lowest standard deviations translate to a higher level of uncertainty, and, ideally, should be covered with a double or triple bet, if possible. The following pseudo-algorithm should describe the overall workflow:               for each match, sorted by std. dev         // "uncertain" matches first         if still can make triple bets             mark top 3 symbols of match        // mark 3 (all) symbols         else if still can make double bets             mark top 2 symbols of match        // mark 2 (highest) symbols         else if can only make single bets      // always does             mark top symbol of match           // mark 1 (highest) symbol      So far so good, but I need to tell the algorithm how much I want to spend. Lets say a single bet costs `1` in whatever currency, the formula to calculate how much a multiple bet costs is:               2^double_bets * 3^triple_bets * cost_per_bet (= 1)      Obviously, the algorithm should try to allocate as much money available as possible into the bet suggestion (it wouldn't make much sense otherwise), and now is where this gets trickier... Lets say I wanna pay a maximum of `4`, listing all possible multiples in PHP (@ IDEOne):               $cost = 1; // cost per single bet     $result = array();     $max_cost = 4; // maximum amount to bet          foreach (range(0, 3) as $double)     {         foreach (range(0, 3) as $triple)         {             if (($double + $triple) <= 3) // game only has 3 matches             {                 $bets = pow(2, $double) * pow(3, $triple); // # of bets                      $result[$bets] = array                 (                     'cost'      => $bets * $cost, // total cost of this bet                     'double'    => $double,                     'triple'    => $triple,                 );                      if ($result[$bets]['cost'] > $max_cost)                 {                     unset($result[$bets]);                 }             }         }     }          ksort($result);      Yields the following output:               Array     (         [1] => Array             (                 [cost] => 1                 [double] => 0                 [triple] => 0             )              [2] => Array             (                 [cost] => 2                 [double] => 1                 [triple] => 0             )              [3] => Array             (                 [cost] => 3                 [double] => 0                 [triple] => 1             )              [4] => Array             (                 [cost] => 4                 [double] => 2                 [triple] => 0             )     )      * * * ## The Problem If I choose to play the maximum amount of money available (`4`) I would have to bet with two doubles, if I use the pseudo-algorithm I described above I would end up with the following bet suggestion:   * `match #1` => `X1`   * `match #2` => `2`   * `match #3` => `12` Which seems sub-optimal when compared to a triple bet that costs `3` and covers more uncertainty:   * `match #1` => `X12`   * `match #2` => `2`   * `match #3` => `2` The above example gains even more relevance if you consider that `match #3` odds could be:               $game['match #3'] = array // stdev = 0.4714     (         '1' => 0.0,           //   0%    home wins         'X' => 0.0,           //   0%    draw         '2' => 1.0,           // 100%    away wins     );      In this case I would be wasting a double for no good reason. Basically, I can only choose the **biggest (possibly stupid) bet** and not the **smartest, biggest bet**. I've been banging my head against the wall for some days now, hoping I get some kind of epiphany but so far I've only been able to come up with two half-[bad-]solutions: * * * ### 1) Draw a "Line" Basically I would say that matches with a stdev lower than a specific value would be triple, matches with a stdev a big higher would be double bets and the rest single bets. The problem with this, of course, is finding out the appropriate specific boundaries - and even if I do find the perfect values for the "smartest" bet, I still don't know if I have enough money to play the suggested bet or if I could make a even bigger (also smart) bet... * * * ### 2) Bruteforce I came up with this idea while writing this question and I know it won't make perfect sense in the context I described but I think I could make it work using somewhat different metrics. Basically, I could make the program suggest bets (# of triple and double bets) for every possible amount of money I could play (from `1` to `4` in my example), applying the pseudo-algorithm I described above and calculating a global ranking value (something like `% of symbols * match stdev` - I know, it doesn't make sense). The bet with the highest ranking (covering uncertainty) would be the suggested bet. The problem with this approach (besides the fact that it doesn't make any sense yet) is that the games my program is going to work with are not limited to 3 matches and the number of double and triple bet combinations for those matches would be substantially higher. * * * I feel like there is an elegant solution, but I just can't grasp it... **Any help cracking this problem is greatly appreciated, thanks.** * * * There seems to be some confusion regarding my problem, I've already addressed this in this question and also in the comments but the misinterpretation still seems to prevail, at least to some. > I need to know **how many triple, double and single bets I will play for a > specific game** (all matches). I already know what symbols I want to play by > looking at each match individually.


title: Two-way ANOVA main effect and interaction query

text: I have a MINITAB printout I was hoping someone could help me with. My question is- from the table below **how do I distinguish between a main effect and an interaction**? Study: Three drug treatments are compared at three sites. The treatments consist of two dosages (low and high dose) of an experimental drug and a control drug. The results were analysed by two-way ANOVA with interaction. The ANOVA is shown below (formatting improved).               Analysis of Variance for response          Source       DF         SS         MS       F      P     site          2     21.301     10.651    0.47  0.656     dose          2    185.404     92.702    4.09  0.108     site*dose     4     90.654     22.663    7.05  0.000         Error        63    202.574      3.215     Total        71    499.933


title: How to test if two factors interact and test if their interaction is important?

text: In a two-way ANOVA, from Kutner's Applied Linear Statistical Models: ![enter image description here](http://i.stack.imgur.com/8teDH.png) What differences are between testing if two factors interact and testing if their interaction is important? I know that testing if two factors interact is the F test via the ratio between the mean interaction sum of squares and the mean error sum of squares. I was wondering how to test if their interaction is important?


title: Bayesian inference with Gaussian distributions

text: This is Problem 4(c), Chapter 2 from  Thrun's Probabilistic Robotics . Note that this is self-study and _not_ homework. Suppose I know my position $x$ to be a normal distribution with density $N(\mu_1=1000,\sigma_1^2=900)$. I query a position sensor, which outputs a measurement $z=1100$; however this sensor is faulty, and its output is a normal distribution with variance $\sigma_2^2=100$. > Given the prior probability and sensor error, how likely is it that the > measurement was the given value,$z$ = $1100$? Given that $p(x) = N(\mu_1,\sigma_1)$ and $p(z|x) = N(\mu_2=1100,\sigma_2)$, the updated position is given by Bayes' rule: $p(x|z) = \frac{N(\mu_1,\sigma_1)N(\mu_2,\sigma_2)}{\int_{\mathbb{R}}N(\mu_1,\sigma_1)N(\mu_2,\sigma_2)dx} $ I found that $p(x|z) = N(\mu=1090,\sigma^2=90)$. But I am unable to take this reasoning further, and am unclear of whether I made any progress. Can you help?


title: Sufficient statistics, MLE and unbiased estimators of uniform type distribution

text: Let $X_1, \dots, X_n$ denote a random sample of size n from the probability distribution with pdf: $$ f_X(x|\theta_1, \theta_2) = \frac{1}{\theta_2 - \theta_1} \ I(x)_{[\theta_1,\theta_2]} \ I(\theta_1)_{(-\infty,\theta_2)} \ I(\theta_2)_{(\theta_1,\infty)}\;.$$ (1) Find a pair of sufficient statistics for $(\theta_1, \theta_2)$. $\bf{My \ thoughts:}$ This wasn't too bad. I got $(X_{(1)}, X_{(n)})$ for this part (2) Find the maximum likelihood estimator $(\hat{\theta_1}, \hat{\theta_2})$ for $(\theta_1, \theta_2)$. $\bf{My \ thoughts:}$ Thinking I need to use monotone functions since it has 2 parameters and the variables are part of the interval. I believe that $\frac{X_{(1)} + X_{(n)}}{2}$ will become one of my estimators. (3) Show that $\frac{X_{(1)} + X_{(n)}}{2}$ is an unbiased estimator for $\frac{\theta_1 + \theta_2}{2}$. $\bf{My \ thoughts:}$ I think I will need to use Cramer-Rao Lower Bound in some form but not quite sure if that is right . (4) Construct an unbiased estimator for $\theta_2 - \theta_1$. $\bf{My \ thoughts:}$ Very stuck on this part, but I think I can use some information from previous parts to help me. Any help is greatly appreciated.


title: Interpreting two-sided, two-sample, Welch T-Test

text: I’m writing some code _(JavaScript)_ to compare benchmark results. I’m using the Welch T-test because the variance and/or sample size between benchmarks is most likely different. The critical value is pulled from a T-distribution table at 95% confidence _(two-sided)_. The Welch formula is pretty straight-forward, but I am fuzzy on interpreting a significant result. I am not sure if the critical value should be divided by 2 or not. Help clearing that up is appreciated. Also should I be rounding the degrees of freedom, `df`, to lookup the critical value or would `Math.ceil` or `Math.floor` be more appropriate?                 /**        * Determines if the benchmark's hertz is higher than another.        * @member Benchmark        * @param {Object} other The benchmark to compare.        * @returns {Number} Returns `1` if higher, `-1` if lower, and `0` if indeterminate.        */       function compare(other) {         // use welch t-test         // http://frank.mtsu.edu/~dkfuller/notes302/welcht.pdf         // http://www.public.iastate.edu/~alicia/stat328/Regression%20inference-part2.pdf         var a = this.stats,             b = other.stats,             pow = Math.pow,             bitA = a.variance / a.size,             bitB = b.variance / b.size,             df = pow(bitA + bitB, 2) / ((pow(bitA, 2) / a.size - 1) + (pow(bitB, 2) / b.size - 1)),             t = (a.mean - b.mean) / Math.sqrt(bitA + bitB),             c = getCriticalValue(Math.round(df));              // check if t-statistic is significant         return Math.abs(t) > c / 2 ? (t > 0 ? 1 : -1) : 0;       }      * * * **Update:** Thanks for all the replies so far! My colleague posted some more info here, in case that affects the advice.


title: What are the practical & interpretation differences between alternatives and logistic regression?

text: A recent question about alternatives to logistic regression in R yielded a variety of answers including randomForest, gbm, rpart, bayesglm, and generalized additive models. What are the practical and interpretation differences between these methods and logistic regression? What assumptions do they make (or not make) relative to logistic regression? Are the suitable for hypothesis testing? Etc.


title: confidence intervals for dependent observations

text: The answer to this question discusses problems associated with calculating P-values for dependent observations. Let's say you have observations from two different groups that are dependent. You consider carrying out a t-test to compare means of the two groups. However, aware that the calculated P-values in a t-test will likely be inaccurate, you calculate means of each group and their 95% confidence intervals. Does calculating confidence intervals avoid problems associated with calculating P-values when observations are dependent?


title: Mean Survival Time: Comparison of Continuous and Discrete Time

text: First, the expected survival time when time is time to event is treated as a continuous r.v. is given as the following: $$\int_0^{\infty} S(t)\, dt$$ The formula for the mean residual survival time $\text{mrl}(u) = E(T -u | T>u)$ is given by $$\frac{1}{S(u)}\int_u^{\infty} S(t)\, dt.$$ Can anyone explain the logic of dividing by $S(u)$? For discrete time, this step does not appear needed (cf related question). Why? Thanks! ADD after Henry's explanation: Lets say I have discrete time (12 months with the following hazards and Survival).               month<-seq(0,12)     h_x<-c(0,0.115706673,0.110186514,0.115769107,0.108296623,0.08908868,0.082548228,0.060146699,0.048112058,0.042197452,0.036919831,0.024691358,0.012787724)     S_x<-c(1,0.8843,0.7869,0.6958,0.6204,0.5651,0.5185,0.4873,0.4639,0.4443,0.4279,0.4173,0.4120)          plot(month,S_x,type="b")      1) Is the mean (truncated) lifetime `sum(S_x[2:12]) #6.3117`? 2) Lets say that we are know that a person survives up to 5th month (start of 5) i.e. they get through 4 complete months and we now want to estimated their mean (truncated) residual life. So, $\gamma$ =7 Is this computed as `sum(S_x[5:12])/S_x[5] #6.358317` since 6.358 < $\gamma$+1 ?


title: When over/under-sampling unbalanced classes, does maximizing accuracy differ from minimizing misclassification costs?

text: First of all, I would like to describe some common layouts that Data Mining books use explaining how to deal with _Unbalanced Datasets_. Usually the main section is named Unbalanced Datasets and they cover these two subsections: Cost-Sensitive Classification and Sampling Techniques. It seems that facing a problem with a _rare_ class you can perform both cost- sensitive classification and sampling. Instead, I think that one should apply cost-sensitive techniques if the rare class is also the target of the classification and a misclassification of a record of that class is costly. On the other hand, sampling techniques, such as over-sampling and under- sampling, are useful if the target of the classification is a good accuracy overall, without focusing on a particular class. This belief comes from the rationale of MetaCost that is a general way to make a classifier cost-sensitive: if one want to make a classifier cost-sensitive in order to penalise a misclassification error of the rare class he should over-sample the other class. Roughly speaking, the classifier tries to adapt to the other class and it becomes specific to the rare class. This is the opposite of over-sampling the rare class, that is the usually suggested way to deal with this problem. Over-sampling of the rare class or the under-sampling of the other class is useful to improve the overall accuracy. Please, it would be great if you confirmed my thoughts. Stated this, the common question facing an unbalanced dataset is: **Should I try to get a dataset that counts as many rare records as other ones?** My answer would be, in case you are looking for accuracy: OK. You can perform it either finding out more rare class examples or deleting some records of the other class. In case you are focusing on the rare class, with a cost-sensitive technique, I would answer: you can only find out more rare class example but you shouldn't delete records of the other class. In the latter case you will not be able to let the classifier adapt to the other class, and the rare class misclassification error could increase. What would you answer?


title: What are the "hot algorithms" for machine learning?

text: This is a naive question from someone starting to learn machine learning. I'm reading these days the book "Machine Learning: An algorithmic perspective" from Marsland. I find it useful as an introductory book, but now I would like to go into advanced algorithms, those that are currently giving the best results. I'm mostly interested in bioinformatics: clustering of biological networks and finding patterns in biological sequences, particularly applied to single nucleotide polymorphism (SNP) analysis. Could you recommend me some reviews or books to read?


title: How to calculate a confidence Interval for a mean of subjects from two sessions?

text: How to calculate confidence interval for this mix of independent(different subjects) and dependent (data of same subject from two sessions) observations? There are 12 subjects and the mean is calculated from these 12 subjects with data from two sessions, so 24 values go into the mean. Any ideas? p.s. R is the language of choice. Thanks!


title: How to read the table to get cumulative relative frequencies?

text: there is a question i can not answer. Suppose the numbers in the following random number table correspond to people arriving for work at a large factory. Let 0,1,and 2 be smokers and 3-9 be nonsmokers. After each ten arrivals, calculate the total relative frequency of smokers and graph the results.               84177  06757  17613  15582  51506  81435  41050  92031  06449  05059       59884  31180  53115  84469  94868  57967  05811  84514  75011  13006       63395  55041  15866  06589  13119  71020  85940  91932  06488  74987       54355  52704  90359  02649  47496  71567  94268  08844  26294  64759       08989  57024  97284  00637  89283  03514  59195  07635  03309  72605       29357  23737  67881  03668  33876  35841  52869  23114  15864  38942        How do I read the table, and get cumulative relative frequency?


title: Survey data analysis in SPSS and in Stata

text: When the survey design is complex, say, a 3-stage stratified cluster sampling design, then we have to use sampling weights. When I use SPSS then I see that the sampling weight it uses is the raw weight. By default it uses Taylor's linearization method for variance approximation. If I want to use relative or normalized weights, then can I get the correct estimate of the variance? The second question is, I have heard that Stata can use normalized weights and some other methods for variance approximation like BRR and bootstrapping method. Can SPSS use these methods? Which one is the best one?


title: Do I apply delta method correctly?

text: I'm no mathematician, so please consider the whole explanation as a subject to criticism and optimization. I'm trying to build a model and calculate an estimate of error variance for the following problem.   Consider a point in 3D, a ray coming out of it towards some plane, and the point of intersection of this ray and the plane. Now the measurements contain the following errors:   1) the plane slope error $\alpha$. I decided to make it isotropic for simplicity and describe it by the angle of aperture of the hypothetical cone, which is drawn by all possible inclinations of the plane's normal vector.   2) the plane shift error $\delta$. This is the parallel shift error (from the fourth member of the plane equation).   3) the direction error $\phi$. I can also be modelled with a cone with an axis coinciding to the error-free ray and with given angular aperture angle. Let's assume that all three random variables are normally distributed and have zero mean (the cones are the same with +- angle aperture, and the shift can go up and down from the error-free state). From some other modelings I also have the estimates of the variances of the three variables. Now I would like to get the the variance of the error in intersection point positioning as a function of relative direction of the ray and position the plane. Of course this error is a random function of the three errors, too, so I have to use the multivariate delta method to estimate the variance of the error given the variances of the input errors. To apply it, I need the gradient vector of the function in the point (0,0,0) as all three variables are zero-mean. To estimate it, I subsequently do the following.   1) Compute the $f(\epsilon,0,0)/\epsilon$ which can be thought of intersection of a slightly inclined plane with no erroneous shift with a error-free ray, which gives one point.   2) Compute the $f(0,\epsilon,0)/\epsilon$ which gives the intersection of un- inclined, but slightly shifted plane with a ray, which again gives one single point.   3) Compute the $f(0,0,\epsilon)/\epsilon$ which is the intersection of a correct plane with a cone which gives an ellipse, from which I can find the longest distance to the "correct" point. Now, having the gradient vector J and the covariance matrix $\Sigma=diag(\sigma(\alpha)^2,\sigma(\delta)^2,\sigma(\phi)^2)$, I can get the error variance as $\sigma(f)^2 = J^T*\Sigma*J$ I do actually get some more or less acceptable values as a result, and even the functional dependence of $f$ on the relative position of the ray and the plane is modeled somehow. But I would like to understand how accurate is this modeling and how viable is my modeling in general. So my questions: 1) if the function $f(\alpha,\delta,\phi)$ is somewhat non-linear, is delta method prone to significant errors? Does it approximate $f$ to be linear with slope $J(0,0,0)$?   2) it seems very unintuitive that in all real cases, when $\alpha,\delta,\phi\neq0$ the error is modeled with an intersection of the cone and a number of shifted and inlined planes, for calculations with delta method I only have to use degenerate cases, as I need the gradient in the vicinity of the $(0,0,0)$ point. 3) does the model and its parameterization look convincing in the least possible order of approximation? I also have some other similar modelings, so I would like to know if my approach is viable. Thanks in advance.


title: Should a mixed effects model be used?

text: This post provides an excellent example of the inner workings of a mixed effects model: http://emhart.github.com/blog/2012/11/16/making-sense-of- random-effects/ In a hypothetical study, I have measured the wing length of a bird species at 10 different locations:               df <- data.frame(wing.length=rnorm(30),                      location=paste('location', rep(1:10, 3)),                       county=paste('county', sort(rep(1:3, 10))))      First 15 rows of dataframe:                  wing.length    location   county     1     29.29024  location 1 county 1     2     30.06387  location 2 county 1     3     27.72127  location 3 county 1     4     29.74502  location 4 county 1     5     29.85506  location 5 county 1     6     30.26669  location 6 county 1     7     30.58748  location 7 county 1     8     30.47608  location 8 county 1     9     31.86882  location 9 county 1     10    28.87578 location 10 county 1     11    30.00726  location 1 county 2     12    30.52488  location 2 county 2     13    30.64339  location 3 county 2     14    29.93695  location 4 county 2     15    27.86217  location 5 county 2      However, unlike the mixed effects model in the linked blog post where each individual is measured 5 times, each location is only measured once. In my hypothetical example, I want to test the hypothesis that wing length varies among counties. County is a fixed effect and location is a random effect. Should a type of mixed effects model be applied, or would a two-way ANOVA be more appropriate? If a mixed effects model is needed, then what type?


title: How do I determine the significance of a difference in weighted averages from a Likert-type survey question?

text: I have multiple survey questions that give respondents 4 Likert-type options for a response to qualify their knowledge/skill in a particular area. I have calculated a weighted average (and divided by the sample size of each group) for the two groups and would like to determine whether there is a significant difference between the two groups. How do I do this?


title: What is the distribution of linear model parameters?

text: I am interested in testing if linear models are statistically different, so I would like to know what I can assume about the distribution of linear parameter models, like the slope for example. More specifically, assume I randomly sample n $(x,y)$ pairs from a random source, and perform a linear regression on these values. What will the distribution of the resulting slope and intercept be? Thanks!


title: Approximate F-test of the intercept with pbkrtest

text: I'm currently trying out the `pbkrtest` package. The test of the fixed effect (`f1`, in the example below) doesn't pose any problem.               ## Not run:     library(lme4)     library(pbkrtest)          # large model:     m1 <- lmer(x ~ f1 + (1|ID) + (0 + f1|ID), data)      When I want to obtain an approximate F-test based on the Kenward-Roger approach for `f1`, I specify this model:               # small model     m2 <- lmer(x ~ 1 + (1|ID) + (0 + f1|ID), data)      and then:               KRmodcomp(m1, m2)      However, what kind of model should I specify if I would like to get an approximate F-test of the intercept?


title: Methods for time-series prediction depending on multiple parameters

text: We have hourly time-series data of the status of a system: number of people present at different train stations. We collected it for a year, and we want to use it to train a model to predict the status of this system in future. We know that several parameters influence the system status, such as: temperature, type of day (holiday/working day/weekend), period of the year, etc. Thus, we collected data for these parameters in parallel to the data of our system status. We are interested in implementing a tool that predicts the status of the system in the next hours based on the previous observations of the system status and on the external parameters that affect it. What methods are most suitable for such a problem?


title: Order of nested random effects in lme4

text: I just fit a model in lme4, and I'm wondering what the heck I fit... I have individuals `id`, and each is measured pass/fail on items that can be described using two factors, `f1` and `f2`. My theory says that `id` and `f2` can interact. So I want to compare these models:               # f2 doesn't matter     m1 <- lmer(pass ~ f1 + (1|id), family="binomial")          # f2 is a fixed effect, doesn't interact with individuals     m2 <- lmer(pass ~ f1 + f2 + (1|id), family="binomial")          # f2 is a nested random effect, interacting with individuals     m3 <- lmer(pass ~ f1 + (1|f2/id), family="binomial")      First, am I right that there's a random effect of each level of f2, and each individual's level of f2 is partially pooled, depending on the number of observations for that individual? Second, does it make sense to fit `f2` as a fixed effect with random individuals nested within each level of `f2`? How does one write that? Third, what's this, with f2 and id reversed in the formula? I typed it by accident, but now I want to understand it.               m3weird <- lmer(pass ~ f1 + (1|id/f2), family="binomial")


title: Bayesian belief network for finding combinatory effects of multiple environmental variables on allele frequency

text: I'm a beginner to Bayesian belief network (BBN). I read a few articles on introduction to BBN. So I know a general idea of BBN. But I'm struggling to construct a graphical network and conditional probability table (CPT) with my data set. A main aim of my research is to find any combinatory (synergistic, additive, antagonistic) effects of multiple environmental variables on allele frequency of some genetic loci. Please go to this link to download a part of my data set. I have 15 environmental variables (each denoted as "var" in the excel sheet) and an allele frequency data (locus1 in the excel sheet) collected at 18 sites. I made a pairwise correlation table of environmental data to see some interactions between the environmental variables but this does not provide information on combinatory effects of them on allele frequency (I think....). I would really appreciate it if you could give me any suggestions on how to make a network graph and CPT with my data set. I use R. So if it's easier to guide me with R scripts, please do so. Thanks so much in advance!