title: On population variable importance

text: Consider we run a random forest on $n$ independent realizations of a random vector $(X_1,X_2,X_3,Y)$ assuming $Y$ is a numerical response variable. Let $f$ be the best theoretical classifier defined by $E[Y\mid X_1, X_2, X_3]=f(X_1, X_2, X_3)$. The _variable importance_ $\hat I_j$ of the $j$-th predictor $X_j$ is defined as follows. For each tree $t$ of the random forest, there's a classifier $\hat f_t$. The mean squared error $MSE_t$ in tree $t$ is the mean squared prediction error in the out-of-bag sample of tree $t$. The $j$-th perturbed mean squared error $MSE_t^{(j)}$ is defined similarly after randomly permuting the values of $j$-th variable in the out-of-bag sample. Finally $\hat I_j$ is defined as the average difference $\overline{MSE}^{(j)} - \overline{MSE}$ over all trees. Let's take $j=1$ to shorten the notations. In Correlation and variable importance in random forests (Gregorutti & al), an expected result is claimed : when $n$ and $ntree$ go to $\infty$, then $\hat I_1$ goes to the _population variable importance_ $$I_1 = E\left[{\bigl(Y-f(X'_1,X_2,X_3)\bigr)}^2\right] - E\left[{\bigl(Y-f(X_1,X_2,X_3)\bigr)}^2\right]$$ where $X'_1$ is a random variable having the same distribution as $X_1$ but independent of all other random variables $X_2,X_3,Y$. Of course some kind of consistency of the random forest classifier $\hat f$ (aggregate of the $\hat f_t$) is required in order that this convergence result holds. ## Question 1: about the proof of the result According to the authors, this convergence result $\hat I_j \to I_j$ is shown in Reinforcement learning trees (Zhu & al). However I only see in this paper some results about the variable importance at a node $A$ of a tree, whose population version is denoted by $VI_A(j)$ in the paper. I don't understand this notion. Using the notations of the paper, a node $A$ is defined by some intervals $(a_j, b_j]$. I'm under the impression that these intervals are treated as deterministic in the paper, whereas (I think) these are the intervals defining the partition of the predictors space, and they are data- dependent in usual random forests. This is the first point I don't understand and I would appreciate clarification. And as a second point I'd like to understand how can one deduce $\hat I_j \to I_j$ from the results of Reinforcement learning trees (Zhu & al). I confess I have not seriously read this paper, maybe $\hat I_j \to I_j$ cannot be directly deduced from the results stated in the paper but can be shown by adapting some proofs. ## Question 2: interpretation of the population importance Correlation and variable importance in random forests (Gregorutti & al) investigates some particular cases such as the additive case $f(X_1,X_2,X_3)=f_1(X_1)+f_2(X_2)+f_3(X_3)$. It is easy to check that $I_j=2Var\bigl(f_j(X_j)\bigr)$ in this case. In the general case, how could we interpret a zero variable importance, a negative variable importance ? How is it related to conditional independence ? ## Question 3: about the `randomForest` R package I have experimentally check that the convergence $\hat I_j \to I_j$ seems to hold with the `rf()` function of the `randomForest` R package. For example:               # simulation      f <- function(x1, x2, x3)  x1 + x2 + x3*(1+x1/5)     set.seed(666)     n <- 20000     sd0 <- 2     x1 <- rnorm(n, sd=sd0)     x2 <- x1*rnorm(n, sd=sqrt(2)*sd0)     x3 <- rnorm(n, sd=2*sd0)     sigma <- 1     y <- f(x1, x2, x3) + rnorm(n, sd=sigma)     # approximate population importances     xx1 <- rnorm(n, sd=sd0)     ( iv1 <- mean((y-f(xx1,x2,x3))^2) - sigma^2 )     ## [1] 12.85     xx2 <- xx1*rnorm(n, sd=sqrt(2)*sd0)     ( iv2 <- mean((y-f(x1,xx2,x3))^2) - sigma^2 )     ## [1] 65.2     xx3 <- rnorm(n, sd=2*sd0)     ( iv3 <- mean((y-f(x1,x2,xx3))^2) - sigma^2 )     ## [1] 36.35     # random forest      ( rf <- randomForest(y ~ x1+x2+x3,  mtry=2, ntree=500, do.trace=100, importance=TRUE) )     ##      |      Out-of-bag   |     ## Tree |      MSE  %Var(y) |     ##  100 |    1.398     2.50 |     ##  200 |    1.373     2.45 |     ##  300 |    1.369     2.45 |     ##  400 |    1.361     2.43 |     ##  500 |    1.362     2.43 |     ##      ## Call:     ##  randomForest(formula = y ~ x1 + x2 + x3, mtry = 2, ntree = 500,      do.trace = 100, importance = TRUE)      ##                Type of random forest: regression     ##                      Number of trees: 500     ## No. of variables tried at each split: 2     ##      ##           Mean of squared residuals: 1.362     ##                     % Var explained: 97.57     importance(rf, scale=FALSE)[,"%IncMSE"] / c(iv1,iv2,iv3)      ##     x1     x2     x3      ## 1.0008 0.9364 0.9703      _But be careful that the convergence can be very inefficient in many situations._ My question is: why the importance measure is named `%IncMSE` ? Since there's the symbol `%` I thought that it could be $100 \times \frac{\overline{MSE}^{(j)} - \overline{MSE}}{\overline{MSE}}$ instead of $\hat I_j=\overline{MSE}^{(j)} - \overline{MSE}$. This is an increase and not a percentage increase, so why this name ?


title: Significant main effect in one-way, but not in two-way, MANOVA

text: I used a one-way MANOVA to study the effect of age groups on the average of height and weight and found that they were significant. Then I used a two-way MANOVA to study the effect of age groups and gender on the averages of height and weight and found that gender had an overall effect but not age groups. Is there any logic in the result obtained using the above two methods? **Update.** The data are of 75 males and 75 females and I randomly made age groups 18-25, 25-30, 30-35, 35-40 and 40+. The number of members in each group varies.


title: Data and model blocks in JAGS

text: I am trying to model a function of observed data in JAGS. For example                  x[i,t]/phi ~ dpois(mu)      where x is observed and phi is a parameter in the model. I can't seem to get this working and am wondering if I may be missing some larger issue in JAGS. I am pretty sure you can't have any transformation of the data in the model statement and need to use a data block but that doesn't seem to be dynamic, thus not allowing for me to fit phi. Any ideas?


title: Supervised or Unsupervised Clustering

text: I have a set of N samples belong to K classes. I am using k-means clustering with Euclidean distance in order to cluster these samples into K clusters. To help the k-means algorithm to group samples of a specific class into one cluster, I initialized the k-means algorithm so that the mean value of each cluster is the mean value of the samples of a specific class. My question is **not** about the usefulness of what I mentioned above or the purpose of it, however, my question is: In machine learning, is this called "supervised clustering", "semi-supervised clustering" or just "normal clustering initialized by means of samples of the real class"? I want to know the correct terminology when we use k-means initialized by real-class mean values.


title: Why were results from PROC MIXED same as PROC GLIMMIX?

text: I am running a simple code comparing insect abundance on a vegetated and unvegetated surfaces. Design is RCBD. 1 factor with 2 level (vegetated and unvegetated). 10 blocks. Distribution of residuals is non normal and variance of the two populations are unequal. I could log transform but wanted to explore using GLIMMIX instead. I obtain the same p values, standard erros and variance using either Proc mixed or proc glimmix. Any idea why? Anything to do with robustness of proc mixed assumptions? I am a lowly environmental student with 2 semesters of stats under my belt...enough to make me confused! Any ideas appreciated.               data abundancevisit3;     input blk trt$ y;     lines;     1 g 5.867     1 b 0     2 g 6.933     2 b 0.444     3 g 0.8     3 b 0     4 g 5.2     4 b 0.667     5 g 56.267     5 b 1.333     6 g 14.933     6 b 0     7 g 54.133     7 b 0.444     8 g 5.026     8 b 0     9 g 4.8     9 b 1.333     10 g 11.733     10 b 0     ;     options nocenter;     proc print data=abundancevisit3;     run;          proc glimmix data=abundance visit3;     CLASS blk trt;     MODEL Y=trt;     random blk;     lsmeans trt;     run;


title: Ellipse region shape from bivariate normal distributed data?

text: In my previous question I needed to help with ellipse region extraction and determine if point lies in that region or not. I ended up with this code:               library(ellipse)     library(mvtnorm)     require(spatstat)          netflow <- read.csv(file="data.csv",head=FALSE,sep=" ")     #add headers     names(netflow)<-c('timestamps','flows','flows_tcp','flows_udp','flows_icmp','flows_other','packe ts','packets_tcp','packets_udp','packets_icmp','packets_other','octets','octets_tcp','octets_udp','octets_icmp','octets_other')     attach(netflow)          #load library     library(sfsmisc)     #plot     plot(packets,flows,type='p',xlim=c(0,500000),ylim=c(0,50000),main="Dependence number of flows on number of packets",xlab="packets",ylab="flows",pch = 16, cex = .3,col="#0000ff22",xaxt="n")     #Complete the x axis     eaxis(1, padj=-0.5, cex.axis=0.8)          pktsFlows=subset(na.omit(netflow),select=c(packets,flows))     head(pktsFlows)     #plot(pktsFlows,pch = 16, cex = .3,col="#0000ff22")          cPktsFlows <- apply(pktsFlows, 2, mean)     elpPktsFlows=ellipse::ellipse(var(pktsFlows),centre=cPktsFlows,level=0.8)          png(file="graph.png")     plot(elpPktsFlows,type='l',xlim=c(0,500000), ylim=c(0,50000))     points(pktsFlows,pch = 19, cex = 0.5,col="#0000FF82")     grid(ny=10,nx=10)     dev.off()          W <- owin(poly=elpPktsFlows)     inside.owin(100000,18000,W)      This produces this graph. ![graph ellipse](https://github.com/matejuh/doschecker_wiki_images/raw/master/linear_regression/density/ellipse.png) Here is the same data with the regression line plotted ![Plot all with linear regression line](https://github.com/matejuh/doschecker_wiki_images/raw/master/linear_regression/density/linRegAll.png). Can you explain me, why the ellipse has this shape? I expected that main axe of ellipse will have the same direction with linear regression line, but it hasn't. Btw. kernel density estimation also points to 100000 althought there are no points... ![kernel density estimation](https://github.com/matejuh/doschecker_wiki_images/raw/master/linear_regression/kernel/kernelPoints.png)


title: Why am I not getting a symmetrical correlation matrix in R?

text: I am trying to obtain the correlation between two matrices of equal dimensions. I have used both column by column correlation and the canonical correlation but I don't get a symmetrical correlation matrix. I get something like the following:                    A    B   C     A  0.20 0.51 0.48     B  0.21 0.56 0.43     c  0.25 0.58 0.45      Any ideas why? This is how I generate the data:               mat1 <- matrix(rnorm(100),10)     mat2 <- matrix(rnorm(100),10)      And this is how I calculate the correlation               correl <- matcor(mat1,mat2)


title: Confidence intervals for 6-point Likert scale

text: I was asked to compare certain ratings between two different groups. The distribution of ratings looks something like this: ` Value Count Percent 1 4 1.32% 2 5 1.65% 3 14 4.62% 4 48 15.84% 5 173 57.10% 6 59 19.47%` 1 is "Very bad", 6 is "Very good". Would it be acceptable to compute confidence (or Bayesian probability) intervals as if the data were normal? If not, what's a better solution. The people who I report to are very non-technical.


title: Effect of including/excluding predictors on beta estimation in linear regression

text: I'm interested in the effect on beta estimation of including/excluding independent variables in linear regression. I've made this data below:               set.seed(50)     predictor1 <- rnorm(10, 3, 1)     predictor2 <- rnorm(10, 6, 1)      I've also simulated a model using these data, with an intercept of 2, a beta of 50 for `predictor1` and a beta of 2 for `predictor2`:               response <- 2 + (50 * predictor1) + (2*predictor2)      A linear regression in `R` correctly calculates the intercept and both the betas:               summary(lm(response ~ predictor1 + predictor2))          Call:     lm(formula = response ~ predictor1 + predictor2)          Residuals:            Min         1Q     Median         3Q        Max      -4.150e-14 -6.833e-15  1.100e-15  9.879e-15  2.757e-14           Coefficients:                  Estimate Std. Error   t value Pr(>|t|)         (Intercept) 2.000e+00  8.194e-14 2.441e+13   <2e-16 ***     predictor1  5.000e+01  7.730e-15 6.468e+15   <2e-16 ***     predictor2  2.000e+00  1.353e-14 1.479e+14   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 2.081e-14 on 7 degrees of freedom     Multiple R-squared:      1,  Adjusted R-squared:      1      F-statistic: 2.095e+31 on 2 and 7 DF,  p-value: < 2.2e-16      However, a linear regression with only `predictor2` in the model calculates a beta of -3.467 for `predictor2`:               summary(lm(response ~ predictor2))          Call:     lm(formula = response ~ predictor2)          Residuals:         Min      1Q  Median      3Q     Max      -75.201 -24.049   5.193  38.423  56.074           Coefficients:                 Estimate Std. Error t value Pr(>|t|)     (Intercept)  171.133    177.588   0.964    0.363     predictor2    -3.467     30.872  -0.112    0.913          Residual standard error: 47.58 on 8 degrees of freedom     Multiple R-squared:  0.001574,  Adjusted R-squared:  -0.1232      F-statistic: 0.01261 on 1 and 8 DF,  p-value: 0.9134      Here is plot of `response` against `predictor1` and `predictor2`, with line of best fit calculated through regression:                library(ggplot2)     ggplot(dat, aes(predictor2 , response )) + geom_point() +        geom_point(aes(predictor1 , response )) +        geom_abline(intercept = 2, slope = 2) +        geom_abline(intercept = 2, slope = 50) +        ylim(0, 300) +        xlim(0,10)      ![enter image description here](http://i.stack.imgur.com/fXEE8.jpg) And here is plot of just `response` against just `predictor2`, against with line of best fit calculated through regression:                ggplot(dat, aes(predictor2 , response )) +            geom_point() +            geom_abline(intercept = 171.133, slope = -3.467)       ![enter image description here](http://i.stack.imgur.com/svEXN.jpg) I have two questions:   1. How did the linear regression arrive at a beta of -3.467 in the second model, when the real beta is 2? Or in other words, how has excluding `predictor1` from the second model caused `predictor2`'s beta to drop by over 5?   2. Is someone able to provide a visual display of how linear regression calculates betas?


title: Detecting online aimbot cheaters in a multiplayer game

text: I am currently working on a anti-cheat mechanism for a multiplayer client called sa-mp. Problem is, recently sa-mp released a new update that fixes latency and hence allows for aimbotters. Aimbotting can be defined as using a computer program to automatically target and shoot all enemies. Problem is, how can aimbotting be detected from the server end? When a player himself is looking at another player's screen, aimbotting becomes obvious. However i require some formula to calculate the probabilities of a player using aimbot so that i can redirect dead players / administrators to spectate the suspected player. While the server cannot give me all the data required to check aimbotters (I did not program the server, i develop scripts to run on it), some relevant available data includes the recent kills by the player, XYZ coordinates of the player's camera and the vectors of the player's camera along with recent bullets fired (And where they hit, whether they hit a player or not e.t.c). I am still in 11th class (10th grade by age i believe), i am yet to study mathematics that can help determine the probability of someone using aimbots. I was hoping somebody could help me out.


title: Can a neural net (MLP) approximate a GLM?

text: I know the activation function (logistic) is linear around 0. Does this then imply a neural net should be able to approximate a GLM given the right bias and weights? I am currently playing in R with a simple model `y~x1`. The GLM performs well out of sample, yet the nnet with 1 hidden variable does poorly. What am I missing?


title: What to do for Dirichlet distribution when elements of X vector may be zero

text: So I want to define a Dirichlet distribution over frequency vectors which are unit vectors whose elements represent the frequency with which different characters occur in a body of text. Trouble is not all the frequency vectors are guaranteed to use every character so some of the elements of the vectors can be zero. This seems really bad when you consider that you will get things like $ 0^{a_i-1} $ in the product which would be a divide by zero error if ${a_i-1}$ is negative. Is there a convention for handling this? I was considering just saying that $x < 0 \implies 0^{x} = 1 $. That way I effectively just drop elements that would be divide by zero errors; im not sure if this is ideal however.


title: Interpretation of $R^2$ in fixed-effects panel regression

text: For cross-section analysis, $R^2$ represents the fraction of the variation that is explained by the model. I understand that a fixed-effects panel regression is designed to optimize for the "between $R^2$". However, what is it's _interpretation_? What is the interpretation of the "within $R^2$" and the "overall $R^2$" (even though they are not optimized for, I assume they still carry some meaning)? Thanks!


title: Is it normal to obtain better (smaller) P values in multivariate analysis compared to bivariate one?

text: If a multivariate design controls for other predictors when calculating the effect of a predictor, shouldn't it give paler P values (less significant ones, or less vivid odds ratios)? I am seeing quite the opposite. Is it **normal** or **possible** or **usual**? Detailed explanation on the case: When I analyze the the correlations between my outcome variable and my five predictors using bivariate correlation coefficients, sporadic few significant P values emerge (all only significant at 0.05 level). However, modeling the same variables within a multivariate logistic regression analysis gives me lots of significant P values, many of which are highly significant (at 0.01 level). I should add that I have modeled independent variables' interactions as well (only their 2-sided interactions). But even if I do not add the interaction terms to the model, still I am getting better results with the multivariate analysis. I should add that none of the variance inflation factors (VIFs) are greater than 2, and I am rather confident multicollinearity is not disrupting my model. So it is interesting to see a multivariate analysis is giving better results. So I wonder is it OK? Or perhaps it is the way it should be (meaning that on the contrary to my belief, a multivariate analysis usually gives better results).


title: Gene expression statistical analysis

text: I need to confirm statistical difference for gene expression data, for two groups of data, standardized with different reference genes:               Sample/reference gene    YLS8/DRH1   14 ref. genes     10C+                   266.5474006   249.6696236     10C-                   11.41636521   9.09824428     20C+                   243.241737    200.5418459     20C-                   5.650170237   5.788838867     30C+                   118.0172933   135.2698666     30C-                   9.867119719   9.964417208      Which way can I do it?


title: Simple statistics for measuring binary time series characteristics

text: Which statistical metrics can be used for characterizing binary time series? For the classical time series there is this paper: > S. Wang and R. Hyndman "Rule induction for forecasting method selection: > meta-learning the characteristics of univariate time series,” > Neurocomputing, 72, 2581–2594, 2009 But for binary time series I do not know any publication.


title: Logistic regression with few observations

text: I suppose regression analysis with few observations has some peculiarities. Here are results for my logistic regression. $R^2_{\rm LR}$ = 0.7, Number of 1 = 6 (4 predicted), Number of 0 = 78 (all predicted) Is it good-fitting? If the number of observation is 10 times higher I would have on doubt but here I'd like to ask for an expert's opinion.


title: Is cosine similarity a classification or a clustering technique?

text: In document classification, is cosine similarity considered a classification or a clustering technique? But you need training data with the cosine similarity for creation of the centroid right?


title: Effect sizes in within-subject design on R

text: I'm doing a meta-analysis but I need your precious help on couple things : 1) First, on R, in the `compute.es` package, I found only effect sizes computation formula for between-subject design. Is there any package/formula to compute effect size for within-subject study? Can we impute the pretest- post-test correlation directly on the same script? As a reminder the variance of d in a within-subject study is (by default, `r=0.5` and does not affect the variance, but it should be interesting to see if there is big difference in the effect size with `r=0.2`, multiplying the variance by `1.6`, or r=0.8 --> V*0.4) : $$V_d=\left\\{\frac{1}{n}+\frac{d}{2n} \right\\}2(1-r)$$ 2) If I have nested (clustered) measures (like, for example, alcohol and cigarette consumption on the same participant), is there a package to do one global effect sizes including those different DV's? 3) For the global effect across studies, I'm using `metacum` in the `meta` package on R. Am I right or is there a more relevant script to do that?


title: Mapping of input and output values for a Multilayer Perceptron

text: I try to understand Multilayer Perceptrons for function regression. As a starting reference I am using the matlab neuronal network toolbox. There I mentioned that the tool maps the input and output values to the range of [-1 1]. Why do they do this? I thought that the inputs need to be decorrelated and free of mean, but they don't do this.


title: Help with regression and GPA

text: I'm currently analyzing some data to establish if school performance for particular subjects makes for good predictors of university performance at first year level across all faculties and degrees. In essence we want to know which specific school subjects are the best predictors of university performance so those can be used for admission purposes alongside admission tests. This has obviously been done a million times before. The problem is that, due to recent curriculum overhauls, I only have access to one cohort's data for now. For various reason, some discussed below, regression the way I was taught might not be the way to go. My first instinct was to do some hierarchical regression for each Faculty by entering the compulsory school subjects needed for admission (for example, math and science for engineering) first and then entering some other subjects of interest second (history and economics for engineering). The outcome variable would be 1st year GPA. This is probably a bad idea though for a number of reasons. 1) GPA, as a construct, will vary from degree to degree within faculties. For example, it would be very misleading to equate the GPA for a Human Resource Management degree and an Actuarial Science degree with one another, though they are both in the same faculty. Theoretically this is a problem but practically as well. For example, despite very similar mean 1st year GPA scores, the mean school math and science mark for Actuarial Science is about 25% higher than the HR Management means. Regression won't like this very much. As far as I can see there's no real way to quantify how much more difficult certain clusters of university subjects are than other clusters making GPA sort of useless as an outcome variable at faculty level. 2) Ideally I'd like to do this at the level of individual degrees where there is some stability in terms of degree difficulty, but if I do this I end up chopping the student numbers severely. With all said and done I have around 5 degrees out of about a 100 yielding sufficient student numbers to run regression. I played around with some power calculations and with 5 or 6 predictors I'll need somewhere in the vicinity of 100+ students. However, due to the selectivity of some of these degrees, school marks tend to cluster at the higher end with very little variance between students and makes for less than ideal predictors. University GPA scores tend to be more dispersed though. I'm hesitant to run regression on degrees with 80 students or less for example, especially given the real risk of collinearity issues with regards to school subjects. There is a pretty robust correlation between math and science marks in the data for instance and all degrees that require one also require the other for admission. Any ideas? Help would be greatly appreciated as I'd hate to construct 100's of regressions with only a small number being usable.


title: weighted regression with dummy variables

text: I am looking for some references on weighted regression when we have only dummy variables for predictors. I would appreciate it. Thanks.


title: condition for a ARMA process to be wide-sense stationary

text: For a ARMA process,   * some (e.g. in Tsay's Financial Time Series) said: it is wide-sense stationary, iff all the roots of its AR characteristic polynomial are greater than 1 in magnitude. This is shown by requiring the variance of the series at any time to be nonnegative and constant.   * some (e.g. in Brockwell's Introduction to Time Series) said: it is wide-sense stationary, iff all the roots of its AR characteristic polynomial are not equal to 1 in magnitude. It is causal, iff all the roots of its AR characteristic polynomial are greater than 1 in magnitude.  Questions   1. According to Brockwell, it seems that a ARMA process can be wide-sense stationary while non-causal, iff all the roots of its AR characteristic polynomial are less than 1 in magnitude. For example, in AR(1) model $X_t = a X_{t-1} + e_t$, the root of its AR characteristic polynomial is less than 1 in magnitude, iff $|a|<1$. I wonder what $Var(X_t)$ in this AR(1) model is then? Is it still well defined (nonnegative)? How does its autocorrelation function behave as the lag increases? (not decay, but stay at 1?)   2. I wonder how you understand the difference in the necessary and sufficient conditions fpr wide-sense stationarity given by Tsay and by Brockwell? Thanks!


title: Time series from experimental measurement

text: I've got a couple of time series obtained from a small number of replicates,where many variables (all of them continuous) are measured at various time points (longitudinal study). I was looking for a statistical test that could answer whether or not the data supports any type of trend (increasing, decreasing or even oscillations), or alternatively, the actual value stays constant and what I see in my individual or average time series is due to measurement error alone. I've had a look and found the Ljung–Box test - Is is appropriate to answer what I want? If not, what do you suggest?


title: How should I construct a prior distribution with a particular kind of count data

text: For context I will first explain the overall problem that I am working on. I am given a catalog of product names and I am also given a large text dataset that may contain mentions of these catalog products. The problem is that the actual product name in the catalog can be very different from how it is mentioned in the text dataset (therefore regex would not produce very good results). Furthermore, there are many product names (or parts of product names) that are ambiguous which means that they can be present in the text dataset but not necessarily refer to the product in the catalog. The evidence that I am using is the measure of how strong the match is between the catalog product name and the mention in the text (this is done with minimum-edit-distance). I want to try to inject some prior information in order to solve some of these problems. The prior information that I have and think would be useful is a measure of popularity of each product in the catalog. My hypothesis is that popular products are more likely to be mentioned and any candidate mention that is also a very popular product should have its posterior probability of P(catalog_entry | mention_text) boosted and vice versa. I am using wikipedia pageviews of the the wikipedia article associated with the product in the catalog as a proxy for its popularity (note that all catalog products are within the same product type/category so the pageviews should be an accurate representation of relative popularity). So my question is what is the best way to tranform the pageview data into probabilities (i.e. what is the best distribution to model the pageview data with). Pageviews are discrete so I suppose I need a discrete distribution where there is a large probability mass at high values. But pageviews are also unbounded on the right. At first I considered a Poisson distribution but it does really fit with what I want because at large enough values the probability mass will start decreasing. I suppose that I could artificially bound the distribution between 0 and the largest pageview that I observe over all products and then use some kind of Gamma distribution. However, a Gamma distribution is not discrete. Can anyone see a solution to this problem?


title: How do I compare individual likert scale items

text: I am looking to compare individual likert scale items. I measured respondent's level of agreement (in a 5 point-scale) for several items (e.g. from 1 to 5 how much is A competent to treat your condition? how much is B? C? and so on..- 11 items in total). I want to analyse if the scores for each statement are significantly different? I am considering the data as ordinal and as such using non-parametric analysis. So far what I did was a Friedman test (which should indicate if there is agreement between k set of ranking) (p<0.001). So indeed the distribution of the several scores (of the several items) are different. But I also wanted to know which items are statistically significanly different. So I did a Kendall correlation for all items.   * Is this the right approach?    * Can anyone let me know what is the best way to compare likert-scale items?


title: Simulating new x's in regression simulation study

text: One of my homework problems is a simulation that compares three estimators (least squares, ridge regression with known parameters, and ridge regression with estimated parameters) for the following model $$Y_i = \beta X_i + \epsilon_i,\quad \epsilon_i\sim N(0,\sigma^2)$$ I am supposed to do 1000 replications with $X_i\sim N(0,2)$. Initially I generated my $X$ vector of data and used the same $X$ vector for each of the 1000 repetitions (so only thing different between repetitions is what random error gets added on). Then I thought that might be wrong and that I should generate new $X$ data between each repetition. What is the correct thing to do? I can provide code if need be, but not really necessary to answer my question.


title: MANOVA, when independent variable is not multi-level and P value significant in MANOVA but not in ANOVA

text: My research has 4 independent variables and 5 dependent variables. My hypotheses are constructed in such a way that 1 independent variable's effect size is checked on 5 dependent variables. MANOVA is an automatically good choice for this. My independent variables are not built around categories and groups e.g. gender and neither it is multi-level like income level. Hence post hoc tests are not relevant in this case. My concerns are as follows,   1. whether MANOVA'S without post hoc tests is a good choice and as effective? Which values should Ii specifically look upto. Wiki's lamba, partial eta square and between subjects?   2. Secondly is there other alternative statistical test for the above mentioned situation. 3.Thirdly, a grave concern is, when I ran regression analysis and ANOVA of one independent variable on one dependent variable, some P values were insignificant, but in MANOVA those P values were significant. What does it signify, is there a concern here? Ii would be thankful if above mentioned queries are answered Sincere regards Salman


title: Estimating distribution function for dependent observations

text: Let $X_1,\dots,X_n$ be identically but not necessary independent distributed with distribution function $F$. I'd like to estimate $F$ efficiently. In case, $X_1,\dots,X_n$ are i.i.d., we can estimate $F$ by $\hat F_n(x)=\frac1n\sum_{i=1}^n1_{\\{X_i\le x\\}}$. $\hat F_n$ is efficient for estimating $F$ since we can verify easily that $var(\hat F_n)$ is the smallest variance among all unbiased estimators. However, since in my problem, I don't have the i.i.d. condition, so what should I do? Any help will be appreciated.


title: How to compute median in an online fashion?

text: > **Possible Duplicate:**   >  What is a good algorithm for estimating the median of a huge read-once data > set? Imagine you have a large, multivariate dataset that resides on disk. Are there any known methods to efficiently compute median with a minimum number of passes through the data ? I've found a candidate for variance/stddev in the name of Welfod/Knutt algorithm, but what about median ? Thanks


title: What sort of feedback button should I implement to measure the difficulty of a task

text: **\---background info---**   My company has built a Search Engine Optimisation product that can suggest monthly prioritised SEO activities (task suggestions), that will have the most benefit to the end user. For my example, let's say we are interested in 10 keywords and 4 task suggestions. The keywords of interest are: “mortgage, mortgage rate, interest rate, mortgage comparison, loan, loan rates, loan interest, loan interest rate, mortgage calculator, loan repayment calculator”. The possible task suggestions are:   1. Build a link to this page with relevant anchor text. For example, find one website and get them to link to your page so that the link contains the words “mortgage calculator”.   2. Produce some content revolving around the keyword: e.g. make a mortgage calculator plugin.   3. Modify an on-page attribute, e.g. change the page title to contain the word “mortgage calculator”.   4. Generate social media interest: e.g. get people to “like” your mortgage calculator on facebook. With 10 keywords and 4 tasks, there are a total of 40 possible task/keyword suggestions. At the moment the tasks are prioritised using a predetermined “difficulty” rating, based on a poll of search engine optimisation experts. We would like to collect data and develop an algorithm which will update the difficulty ratings. There are two reasons for this:   1. Certain tasks might be more difficult for some clients. E.g. if a client does not have someone in charge of social media, it might be hard for them to implement a social media campaign. Therefore, we would like our difficult ratings to reflect our client's specific circumstances   2. We would like to give new clients better info on how difficult a task is based on our current clients' experiences.  **\---end background info---** I have two questions, one related to data collection, and the second related to how the data collected should be analysed. **Question 1:** We would like to get user feedback on task difficulty. **What sort of button should I implement to poll users** so that we get useful data which can be used to update our difficulty ratings? Should this button be:   1. Something different to the following options?   2. Yes/no response to the question “Was this task too difficult”.   3. Please rate the difficulty of this task on a scale of 1 to 5.   4. Is this task more difficult than the task above/below it?   5. We already have an option for the user to reject the task. Perhaps we could derive information about difficulty from this? **Question 2:** **What sort of analysis methods can we use to update the difficulty rating?** Does the analysis method differ depending on the way we collect data?


title: Performing linear regression analysis on two data series with different sample spacing

text: I have a record of one climate variable with a data point every year, and another one which has sample spacing that varies between 1.3 and 75.2 years. I even have a few ages in that series for which I have two data points. If I want to figure out if these two records are significantly correlated, should I then a) resample the low-resolution record to match the high-resolution one, and then run the regression; b) resample the high-resolution record to match the low-resolution one, and then run the regression; c) do a regression on the two records without any resampling (I don’t know how!) or d) not use regression in the first place (but what then?). Do notice I struggle with statistics, and prefer to limit my data processing to excel (be it with “data analysis add-in”). Thanks a lot.


title: Correct model to use for count data

text: I've been working on an assignment for a while as a part of my final qualification thesis. I'll try to describe my problem in a nutshell. There is a small cable TV channell which sells some types of goods liks blenders, mops, etc. 24/7 they are broadcasting advertisments of the different products (the range is constant, 8 different products) of different duration (5, 10, 15, 20, 30 minutes). The fact of selling is assigned to this particular advertisment if the customer called after its beggining and before the next ad of this product. My task is to construct prediction function for the following broadcast optimization. I need to predict the number of orders for the given advertisment. So, after countless attempts, I chose negative binomial regression to model this process. I have several predictors including interval before the start of the next ad of this product, set of dummy variables indicating products and durations, time of the day, number of this product ads (3 hours before this ad and 3 hours after) and so on. Is there any way to use more suitable model, for example, with fixed or random effects? As I understood, these methods were developed for the panel data but my set doesn't seem to be panel. Thanks for your attention, any ideas would be appreciated. P.S. There is a small sample of data, JIC https://drive.google.com/file/d/0B34-IapABY3mZllqOC1ZeWhWTkU/edit?usp=sharing


title: Sum of Normally Distributed Random Variables

text: I am aware that the sum of two (or more) normally distributed Random Variables is **not** necessarily also normal. I do have questions regarding the special case where I have the following addition: For c $\in \mathbb R$, is c + N(0,a), where I add a normal Random Variable with mean zero and variance a to c, is then the sum normally distributed? Also, does this addition correspond to drawing a random variable X too from a Normal distribution with mean and then add it to N(0,a)? Hope the question is clear, thanks


title: How to run a power analysis for a repeated measures (all within-subjects) with combined contrasts?

text: I am looking for a little help on determining the correct way to run a **power analysis** on a **within-subjects repeated measures** design (i.e., no between subjects variable) with **combined contrasts**. Originally, I was under the impression that I could run a power analysis on a repeated measures design with planned contrasts and not worry about the fact that we are combining the variables because everything is within subject. However, I have since realized that this is not the case because, for a three factor Repeated (Within) design, there will be 7 error terms owing to combination of main effects and interactions. Depending on how you want to form contrasts (and power or effect sizes) we pool different SS error and their associated df terms... A SS Error(A) B SS Error(B) C SS Error(C) AB Error(AB) and so on... For example, we are looking at how different apps impact quality of life measures. One of the planned contrasts is to compare manufacturers (2). However, we would also like to collapse across manufacturers (i.e., combining the planned contrasts) to look at app use versus no app (baseline). This website has been so helpful for me in my research, and I have been struggling for some weeks now to figure this out/get a straight answer out of some wonderful statisticians. Perhaps I am asking the wrong question, but hopefully someone here can help me. Any advice would be much appreciated. Thank you in advance!


title: Is the information matrix equality valid for the Poisson distribution?

text: As far as I know it should, since the support of the Poisson is independent of its lambda parameter. The negative of the expected Hessian equals $\frac{n}{\hat{\lambda}}$, where $\hat{\lambda}$ is the MLE estimator. Information matrix equality therefore implies $E((\frac{\delta logL}{\delta \lambda})^2)|_{\lambda=\hat{\lambda}} = \frac{n}{\hat{\lambda}}.$ But no matter how I turn it, the left hand side of the latter expression always equals zero. This is obviously so because the MLE estimator lets the score disappear (by definition). What am I doing wrong?


title: Assumptions of two way anova

text: Please tell me what the actual assumptions of a two way anova are. I read somewhere that this being similar to multiple regression, the only assumption is that the residuals are to be normally distributed with equal variance. But then a web search gave the following assumption "Two-way anova, like all anovas, assumes that the observations within each cell are normally distributed and have equal variances". I am confused. Please help out.


title: Gibbs sample from AR(1) of exogenous input

text: I am trying to fit a model where there is a sequence of exogenous "shocks", $X_1, X_2, ..., X_T$, and a AR(1) of these shocks explain $Y_1, Y_2, ..., Y_T$. Specifically, Data (known): $X_1, X_2, ..., X_T$ $Y_1, Y_2, ..., Y_T$ Priors (known): $\phi = 2 \phi^* - 1$, where $\phi^* \sim Beta(\alpha, \beta)$, "persistence" parameter, between $(-1, 1)$ $\beta \sim N(\mu, v)$, "shock strength" parameter Model: $E_t = \phi E_{t-1} + \beta X_t$ $Y_t = E_t + \epsilon_t$, $\epsilon_t \sim N(0, 1)$ Unknown parameters: $\phi$, $\beta$, $E_0$ My question is, what is the easiest way to sample $\beta$ and $\phi$? It seems that if I am allowed to use $E_{t-1}$ when sampling from $\phi$, it becomes a bit easier, as: $Y_t = \phi E_{t-1} + \beta X_t + \epsilon$ $\implies \phi E_{t-1} = Y_t - \beta X_t - \epsilon$ $\implies \phi \sim N(\frac{Y_t - \beta X_t}{E_{t-1}}, \frac{1}{E_{t-1}^2})$ However, since $E_t$ is known deterministically from $\phi$, it seems that if I am conditioning on it I must choose a deterministic value for $\phi$. So I have two questions: How do I best fit this model using MCMC? How do I condition on a deterministic auxiliary variable within a Gibbs sampling step?


title: Pattern mining on a small data set

text: I have a small data set 30 features/predictors and 30 observations. My target variable is Oil production and my predictors are well & reservoir properties (depth, trajectory, temperature, pressure etc) I want to identify the manner in which these features impact the target variable via identifying patterns that may exist in my dataset. The relationships may be non-linear and many features are correlated to each other. Can you please suggest me some great ways to approach this problem? Any help will be greatly appreciated! Thanks.


title: ANOVA to find subjects that responded significantly different between treatments

text: I have a dataset composed of several subjects exposed to two different treatments (treat) in three distinct times (test), and I am interested in highlighting which subjects have significantly different response between the two treatments. Example:                  subject treat test    result     1        A     X    1 2.2740321     2        A     X    2 1.5900599     3        A     X    3 0.8478363     4        A     Y    4 1.5506098     5        A     Y    5 2.2662206     6        A     Y    6 1.2903771     7        B     X    1 1.1469458     8        B     X    2 1.5929735     9        B     X    3 1.0960918     10       B     Y    4 1.4392964     11       B     Y    5 1.2450169     12       B     Y    6 1.7282016      Any ideas? Note that I am not interested in seeing which treatment (group) is more effective.


title: MV GARCH DCC - t-test with high t-statistics

text: I use multivariate GARCH DCC to compute volatilities of 3 (stationary, supposedly) time series. Then I test the differences between these three volatility series with t-test and get very high t-statistics. Could anyone help explaining me this? Is there any adjustment I need to use for the t-test?


title: What does the notation $(\textbf{X} \perp \textbf{Y} , \textbf{W}\mid \textbf{Z})$ mean?

text: I was reading Koller's and Friedman's Probabilistic Graphical Models book and became confused about some of its notation because of a set of notes that either contradict it or express it differently. Whichever it is, I would like to clarify. The notation is: $(\textbf{X} \perp \textbf{Y} , \textbf{W}\mid \textbf{Z})$ The reason I got confused on what the notation means is because of the following slides online: http://www.stats.ox.ac.uk/~steffen/seminars/waldmarkov.pdf Particularly what confuses me is the way these notes express the decomposition rule on page 8 and how koller's book expresses this same decomposition rule a little differently. If I understand the two notations, I should understand the decomposition rule under both notations. Koller's book expresses the decomposition rule as follow: $$(\textbf{X} \perp \textbf{Y} , \textbf{W}\mid \textbf{Z}) \implies (\textbf{X} \perp \textbf{Y}\mid \textbf{Z}, \textbf{W})$$ Which I would express intuitively as; the set of random variables $\textbf{X}$ are independent of the set of random variables $\textbf{Y}$ and $\textbf{W}$ (i.e. interpreting the comma in between $\textbf{Y}$ and $\textbf{X}$ as an "and") given the set of random variables $\textbf{Z}$. However, the slides seem to interpret this comma separation as an Union which confuses me (because I thought we were doing "ands"): $$(A \perp_{\sigma} (B \cup D) \mid C) \implies (A \perp_{\sigma} B \mid C) \text{ and } (A \perp_{\sigma} D \mid C) $$ It confuses me because they seem to use union for the commas were I used "ands". Furthermore, I can't see why the two statements are equivalent decomposition rules (which might explain why I've had a tough time proving it because I don't think I understand the statement!). Does someone know the difference in the notation and how they equivalently express the decomposition rule? * * * As a reference this is more background on the notation in Koller's Book: Let capital non-bold stand for random variables say $X$ is a r.v. Let little non bold stand for the assignment to a random variable say $(X = x)$. Also, let me define captial bold letters as sets of random variables. For example $\textbf{X}, \textbf{Y}, \textbf{Z}$ are three sets of random variables. Let small bold letters denote assigments to these sets $\textbf{x}, \textbf{y}, \textbf{z}$ i.e. it denotes assigments of values to the variables in these sets. Let $Val(\textbf{X})$ be the values that the set of random variables can take. On page 21 of Koller's book it explicitly says that Intersections are comma separated, so its fair to assume comma's means "intersection". This can be appreciated on page 21 where it says: > Rather than write $P((X=x) \cap (Y=y))$, we write $P(X=x, Y = y)$ or just > $P(x,y)$. I know this is only on the context of explicit probability distributions, however, if we note how it starts to prove the decomposition property, we can infer that koller probably mean a comma as a intersection in the confusing notation $(\textbf{X} \perp \textbf{Y} , \textbf{W}\mid \textbf{Z})$ too. She starts like this: > For example, to prove Decomposition assume that $(X \perp Y , W \mid Z)$ > holds. Then, from the definition of conditional independence, we have that > $P(X,Y,W|Z) = P(X|Z)P(Y,W|Z)$. Now using... proof continues... The equation: $$P(X,Y,W|Z) = P(X|Z)P(Y,W|Z)$$ Makes me suspect that the slides are either wrong or sloppy in notation. Or both are sloppy. I am not sure.


title: Performing CFA & Structural Models

text: I have developed a conceptual framework based on literature. The conceptual framework comprises antecedents and consequences. The antecedents factor 1, factor 2 and factor 3 affects factor 1 and consequences factor 1 affects factor 4 and factor 5. Since my respondents includes organizations, my sample is around 160. There is slight non-normality in the data. Can I perform two CFAs two structural model - or can I perform CFA and structural model for antecedents and then proceed to CFA and structural model for consequences/outcomes? (The idea to run separate CFA and SM is because relatively small sample and slight non-normality.) Suggestions are appreciated.


title: Retrospective Case Study

text: I am carrying out some research and am a novice to statistics, having done a 2 week module primarily focusing on parametric data and experimental design. However, for my research I am doing a retrospective case study, on the hearing and speech outcomes of 3 groups of patients with 3 different ear conditions following a surgical intervention to improve their hearing. Ear condition being the independent variable, speech and hearing outcomes, the dependant variable. The sample size is based on available data not a power calculation and each group has the following numbers, each group represents a different ear condition: 1) n =10, 2) n= 5 and 3) n=2. I have been advised to do inferential stats, so have done the appropriate tests, but I wonder if inferential stats are geared toward an experimental design? The outcomes are also affected by many other demographic factors between patients such as the age they had the surgery, underlying cognitive ability and so forth- these were not controlled for owing to the type of study. I also wonder as no power was calculated before (normally it is done following a pilot) but would you do a post hoc power calculation in this type of study? I have also done descriptive stats, but I am quite uncertain about the use of inferential stats (Mann Whitney and t-tests) were appropriate to the data...but I am unsure what inferences can be made on this type of study design? Any thoughts would be greatly appreciated.


title: Comparing term-frequency distributions with unequal sample sizes?

text: **Background** I have several datasets of word frequencies where some datasets have much more data than others: from 3000 samples to 20000 samples. I also have large reference corpora with millions of samples. **Question** Using term probabilities, I want to compute the pairwise distance between the datasets and the reference corpora to be able to say "Of all corpora, Dataset $A$ is closest to Corpus $S$ while Dataset $B$ is closest to Corpus $P$". However, I want a metric that will not be overly influenced by the sample sizes of the datasets. For example, if $A$ has 50000 samples, $B$ has 3000, and reference corpora have millions of samples, $A$ might often be closer to the reference corpora than $B$ simply because it has more samples (assuming $A$ and $B$ are similar but not drawn from the same distribution)? I believe that with a naive application of K-L divergence, for example, $A$ would be favored over $B$ but please correct me if this is not the case.   1. What is the best metric for measuring the distance between the term probabilities of a sample and a large corpus?   2. How does unequal sample size affect the metric?


title: Python packages for working with Gaussian mixture models (GMMs)

text: There seem to be several options available for working with Gaussian Mixture Models (GMMs) in Python. At first glance there are at least:   * PyMix - http://www.pymix.org/pymix/index.php Tools for mixture modeling   * PyEM - http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/em/ which is part of the Scipy toolbox and seems to focus on GMMs **Update: Now known assklearn.mixture** .   * PyPR - http://pypr.sourceforge.net/ pattern recognition and related tools including GMMs ... and perhaps others. They all seem to provide the most basic needs for GMMs, including creating and sampling, parameter estimation, clustering, etc. What's the difference between them, and how should one go about determining which is best suited for a particular need? Ref: http://www.scipy.org/Topical_Software


title: Sequential Proportions Test

text: I have been trying to find some R code but haven't had any luck so far. I am able to perform a sequential t-test on some data that I was working on previously but now my data has a combination of real (numeric) and categorical variables in it. I would like to be able to write some code that will take one variable at a time and perform a two proportions test on that column where the samples are defined by a breakpoint that moves sequentially through the observations in that column. Once all the p-values are stored in a matrix the code would move on to the next variable and perform the same testing procedure. At the end of it all I would be left with a matrix of p-values and all that left for me to do is interpret the results. The reason I ask this because I would like to know at what point (i.e. where the breakpoint between each of the samples is) do each of the variables in my sample become significant. I could do it manually but I have over 10 000 observations with 30 variables.


title: Can anyone suggest penalized maximum likelihood method package for OLS type regressions in Stata?

text: I need to use penalized MLE for my work on relationship between distance and its implication for economic outcomes. I have read many theoretical papers but I could find no Stata program/packages to do the analysis. Any help on how to use the penalized maximum likelihood method for OLS type regressions in Stata would be appreciated.


title: is there something like Mann-Whitney U test that can control for a continous variable?

text: I want to compare mean inhibitory concentration (continous variable) in patients with outcome A versus outcome B. I used Mann-Whitney U test. However I would like to control those results for the parasite density, which is also a continous variable. Van Elteren test only allows to correct for strata. Is there any similar test that allows to correct for a continous variable, or does this not make any sense?


title: How to properly group data for analysis

text: I am wondering what the most common way to handle data is wherein a participant belongs to multiple groups. For example, if you're wanting to run linear regression and look at how specific diet-type (IV/predictor) effects a specific health outcome (DV/response variable) and you have a few participants that are on a combination of 2 diet types (lets say "Mediterranean" AND "Vegetarian") and you actually want to know how Mediterranean outcomes are different than Vegetarian, how do you handle those participants? I know you could make another group that's a combination diet type, but not if you have too small of an N for that. Any input is appreciated!


title: Heckman's two stage: multiple selection models

text: If a sample is likely to be self-selected on multiple selection criteria, does it make sense to include Inverse Mill's ratios for multiple selection models in the same second stage OLS model?


title: Stats: Relationship between Alpha and Beta

text: My question has to do with the relationship between alpha and beta and their definitions in statistics. alpha = type I error rate = significance level under consideration that the NULL hypothesis is correct Beta = type II error rate If alpha is lowered (specificity increases as alpha = 1- specificity), beta increases (sensitivity/power decreases as beta = 1 - sensitivity/power) **How does a change in alpha affects beta ?** Is there a linear relationship or not ? Do the ratio alpha/beta always the same, in other words the ratio specificity/sensitivity is always the same ? If yes, it means that by using a bonferroni correction we're just shifting to lower sensitivity and higher specificity but we're not changing sensitivity/specificity ratio. Is that correct to say so ? Update (Case-specific question): For a given experimental design, we run 5 Linear Models on the data. We have a True Positive Rate (sensitvity/power) at 0.8 and a True Negative Rate (specificity) at 0.7. (Let's imagine we know what should be positive and what should not.). If we now correct the significance level using Bonferroni to 0.05 / 5 = 0.01. Can we numerically estimate the resulting True Positive Rate (sensitivity/power) and True Negiative Rate (Specificity) ? Thanks a lot for your help.


title: What does a squared coefficient of variation signify?

text: The coefficient of variation is defined as the ratio of S.D to mean, a dimensionless number which on multiplication by 100 gives the percentage of dispersion with mean. But when we square this coefficient of variation, what does it signify? (I have to understand a graph having squared coefficient of variation on y axis and mean on x axis for Poisson distributed data). Thanks


title: Bayesian $\chi^2$ discrepancy

text: We define the $\chi^2$ discrepancy as $$T(y,\theta)=\sum_i\frac{(y_i-E(y_i|\theta))^2}{\mathrm{var}(y_i|\theta)}$$ (definition from Gelman). How do we apply this if $\theta$ is a random variable, and we have samples of it(but not the distribution)?


title: What is the fiducial argument and why has it not been accepted?

text: One of the late contributions of R.A. Fisher was fiducial intervals and fiducial principled arguments. This approach however is nowhere near as popular as frequentist or Bayesian principled arguments. What is the fiducial argument and why has is not been accepted?


title: How to prove independentness of marginal/conditional (?) posterior distributions?

text: This is a question about exercises 4.2 and 4.3 of Jim Albert’s “Bayesian Comptutation With R” (p. 82). Note that while this might be homework, in my case it is not. We are to prove that, given two indepentent normal distributions `N(µ₁,σ²₁)`, `N(µ₂,σ²₂)` and the prior `g(µ₁,σ²₁,µ₂,σ²₂)∝1/(σ²₁σ²₂)` the resulting posterior distribution is such that the vectors `(µ₁,σ²₁)`, `(µ₂,σ²₂)` are independant. My understanding is the following. Since the distributions are saind to be independent, the joint posterior distribution equals the distribution of `N(µ₁,σ²₁)` times `N(µ₂,σ²₂)` times the prior distribution `1/(σ²₁σ²₂)`. The resulting posterior distribution is: $$g(µ_1,σ_1,µ_2,σ_2|y_1,y_2) ∝ \frac{1}{{σ^{2}_{1}}^{(\frac{N_1}{2}+1)}}exp(-\frac{1}{2σ_1^2}(\sum_{i=1}^{N_1}(y_{1i}-\bar {y_1})^2+n_1(µ_1-\bar{y_1})^2)) * \frac{1}{{σ^{2}_{2}}^{(\frac{N_2}{2}+1)}}exp(-\frac{1}{2σ_2^2}(\sum_{i=1}^{N_2}(y_{12}-\bar {y_2})^2+n_2(µ_2-\bar{y_2})^2))$$ Now, two vectors are said to be independent if the joint distribution equals the product of the marignal distributions. But that is how I designed the posterior density. But the answer „The vectors `(µ₁,σ²₁)`, `(µ₂,σ²₂)` are independent because otherwise my posterior distribution would be wrong“ seems a bit sketchy to me. So there might be another way of proving independence of vectors. The marginal distribution of a continuous variable can be obtained by integrating over the joint distribution [wikipedia]. This is described for the two-variable case. So the marginal distribution of `(µ₁,σ²₁)` would be $$∫_{µ₂,σ²₂}g(µ_1,σ_1,µ_2,σ_2|y_1,y_2)dµ₂σ²₂$$ But I do not know how to compute that integral or how to apply the formula of the two-variable case to my multivariate case, further complicated by the fact that my `g` formula is conditional on `y₁` and `y₂`. So my question is: **Is integration the correct approach to answer that question?** I rather doubt that, because integration as a means to compute marginal distribution has not been introduced by Jim Albert in the book. Instead, he describes a family of _dependent_ priors, that is, Howard’s prior. They can be used in the case of analysis of proportions, not of mean and variance of normal distributions (or so I believe). Howard’s prior may still be interesting for this particular question as we might investigate what mathematical property of Howard’s prior causes dependence. Howard suggests transforming the proportions to logits $$\theta_1 = log \frac{p_1}{1-p_1},\theta_2 = log \frac{p_2}{1-p_2}$$ which enter the prior density in a multiplicatve term $$e^{-(1/2)u^2},u=\frac{1}{σ}(\theta_1-\theta_2)$$ So, could I say that because the difference (or a function thereof) of the random variables enters the equations (first the prior, then the posterior distribution), they are not independent anymore? And how would that help me in answering the original question? Related questions. Is my theory on independence correct? What is the difference between the marginal and the conditional distribution? Am I right in assuming that the marignal distribution is the mean of the conditional distribution, that is, conditioned on all other variables?


title: How to fit data with nonlinear partial least squares in R?

text: I am looking for a way to do nonlinear partial least squares in R or matlab. I thought kernel pls was a way to do it but it is not directly related to nonlinear pls. Do I have to calculate my own nonlinear kernel matrix and pass it to kernel pls? Any hint would be appreciated. Update: I don't think there is any direct way to do nonlinear pls in R or matlab. But I found a package kernlab that can be useful when you deal with kernels for nonlinearity. I think this can be a good start point.


title: Clarification on using spec.pgram in R

text: Does anyone understand how the argument **"span"** is used in **spec.pgram** in **R**? From the help menu, it provides the following examples:               spectrum(ldeaths)     spectrum(ldeaths, spans = c(3,5)) gives bandwidth= 0.241     spectrum(ldeaths, spans = c(5,7)) gives bandwidth= 0.363     spectrum(ldeaths, spans = 3)      gives bandwidth= 0.127        1. What do **c(3,5)** and **c(3,7)** or **spans=3** represent?    2. How the bandwidths are calculated? Thanks!


title: Simulate data for stratified Cox model

text: I need to generate simulated data for a conditional logistic regression using `R`. The code I wrote appears correct, but something is wrong. After generating my data, betas, and labels, I use `glm()` to attempt to estimate the betas. If everything is correct, the estimated betas from `glm` should be very close to the true betas. However, `glm()` fails to converge. Below is my code. Can anyone point out what I'm doing wrong? Thanks!               library(survival)          # Initialize parameters     k <- 5     n <- 10     g <- 20          p <- rep(0, n*g)     y <- rep(0, n*g)          # Create the group index     GI <- vector(mode="numeric")     for(group in 1:g){         GI <- c(GI, rep(group,n))     }          # generate the data     betas <- matrix(rnorm(k), ncol=1)     x <- matrix(rnorm(n*g*k), ncol=k)          # calculate the true labels     for(group in 1:g){         p[GI==group] <- exp(x[GI==group,] %*% betas) / sum( exp(x[GI==group,] %*% betas) )         y[n*(group-1)+which.max(p[GI==group])] <- 1     }          cl <- clogit(y ~x + strata(GI))


title: Order statistics of independent NOT identically distributed random variables

text: Can I find the p.d.f of the order statistics (min for example) from a set of independent, but not identically distributed random variables? (the analytical p.d.f. of the other variables is at hand)


title: Parameter space exploration

text: _I do realise this question is quite specific and practical, but I seek for some general help which helps me progress further in my analysis._ Let $y(\boldsymbol{x})\in\mathbb{R}$ be the function I'd like to regress, whereas $\boldsymbol{x}\in\mathbb{R}^6$. My design matrix $\boldsymbol{X}\in\mathbb{R}^{11\times6}$ (therefore, at the moment, $y(\boldsymbol{X})\in\mathbb{R}^{11}$, but I can run further experiments, if needed).   My aim is maximising $y(\boldsymbol{x})$ over $\boldsymbol{x}$, $x_{2k+1}=\\{5,7,9\\}$, $x_{2k}=\\{2,4\\}$. If I would run a _dense_ experimental analysis, I would to test $3^3\times2^3=216$ cases, and each of those are very time consuming. The collected data is the following (yes, I know, they do not respect the constrain given above, which it allows me to run a reasonably lower total number of experiments if the _dense experimental analysis_ will be considered as only feasible solution):               X = [          4     2     3     2     3     2          7     2     5     2     5     4         10     2     7     4     5     2         10     2     9     2     8     4         10     4     6     2     4     2          7     2     5     4     5     2          7     4     5     2     5     2          7     2     5     4     5     4          7     4     5     2     5     4          7     4     5     4     5     2          7     4     5     4     5     4          ];          y = [          64.7          68.0          69.4          68.5          68.5          71.9          71.6          71.7          71.0          71.6          68.6          ];      I did scatter plot $(x_i, y(x_i)),i=\\{1,2,3,4,5,6\\}$, but the output is not that informative.               name = {          'ks1'          'ss1'          'ks2'          'ss2'          'ks3'          'ss3'          };          [m, n] = size(X); % # of experiments, # of parameters = 3 conv kernel size + 3 subsamp          for j = 1:2:n         subplot(2,3,(j+1)/2)         h = plot(X(:,j),y,'x','MarkerSize',10,'LineWidth',2);         xlabel(name{j})         ylabel('Testing accuracy [%]')         axis([2 10 60 75])         grid on     end          for j = 2:2:n         subplot(2,3,3+j/2)         h = plot(X(:,j),y,'x','MarkerSize',10,'LineWidth',2);         xlabel(name{j})         ylabel('Testing accuracy [%]')         axis([1 5 60 75])         grid on     end      ![Scatter plots](http://i.stack.imgur.com/R7WP6.png) So, my question is, **how can I understand** (and/or visualise) **what is going on in this multivariate problem?**   To my mind a PCA keeps coming up as possible useful tool, even though I know it is used for visualising main orientations rather than helping me finding the maxima of my scalar fields.. **Update**   I am currently running the 216 cases batch experiment, limited to 30 epochs training. Still, I would like to _visualise_ my results. Actually, I do have a testing (and training) accuracy function per case, whereas above I just picked its maximum. This is an overview of the testing accuracy functions ![Testing accuracies](http://i.stack.imgur.com/0mYvb.png) * * * What I am trying to accomplish is determining the best configuration for a convolutional artificial neural network. The parameters I am considering here are _kernel size_ (`ks`) and _subsamplig_ (`ss`) factor for multiple layers (`1`, `2` and `3`), whereas $y(\cdot)$ is the testing (or cross-validation) accuracy.


title: Why is Sampling Importance Resampling (SIR) better than Importance Sampling (IS)?

text: From what I understand, SIR is a mechanism for sampling from a distribution $p$ that works as follows:   1. Approximate a target distribution $p$ using an importance sample $S$ from a proposal distribution $q$   2. Draw a **small** sample $S_\text{small}$ from $S$ with replacement according to $p$ The result, $S_\text{small}$, is a sample from $p$.   1. Why use a **small subset** and not the **full** sample $S$ in 1?   2. Why is this a variance reduction technique, and why does it matter to reduce the variance?


title: Why use bayesglm?

text: My overall question is: why use `bayesglm` instead of other classification methods? Note:   1. I'm only interested in prediction.   2. I have a decent amount of data (~ 100,000 obs).  I feel like the sample size is large enough the parameters of a regular logistic regression are going to be normally distributed (CLT). What would I gain by specifying priors? My hunch is that it will only matter for a small dataset, but I don't have any theoretical or applied evidence.


title: Finding class probability for a set based on subsample classification

text: I have 2 sets of X typed objects:               A = {x1, x2, x3 ... xn}     B = {x1, x2, x3 ... xm}      The goal is to determine how similar those two sets are based on the output of a probabilistic random forest. This random forest takes as input two of these objects (one from each set) and then returns the probability that they belong to the same class. So far so good. The classification process is expensive so instead of computing every possible combination of elements between **A** and **B** , I randomly select **k** pairs of elements to be classified, so now I have a vector of **k** elements, each representing a probability of that pair of elements belonging to the same class.               {P1, P2, P3 ... Pk}      **How can I find the right _k_ so it is a statistically representative set of samples, but small enough to reduce computational cost?** Based on that, I need to determine the probability of set **A** and **B** belonging to the same class. **How can I compute the probability for the two sets based on the set of samples taken?**


title: Reproduce boxplot with only data points at given quantiles

text: It is trivial to create a boxplot in R with a full dataset. However, with limited access to the whole dataset, I just have 5 data point at min, 25%, 50% ,75%, and max. So is there any easy way to reproduce the boxplot with only these 5 values?


title: Ratio data with categorical events

text: I have a question about treating qualitative (categorical) events within otherwise quantitative (ratio scale) data. Without going into too much detain, the experiment is on throwing. I measure the distance between the target and the result of the throw for each attempt. Some attempts are failed because the participant hits an obstacle on the way to the target. In the end the dependent variable may look something like this: [.2 .1 .3 .2 -1 .2 .2 .8 -1] where the values may range from 0 to anything, but -1 represent the obstacle hits. What are some ways to deal with such categorical data within quantitative variable? I cannot just throw them out because they are meaningful from the performance perspective. Substituting with some high error distance seems arbitrary. Is there anything else? In the end I'm interested in how people get better at throwing as a result of practice.


title: Testing if alcohol consumption and smoking are independent

text: The question asks to test if smoking status and level of alcohol consumption are independent using the usual five-step procure at alpha $=0.05$: ![enter image description here](http://i.stack.imgur.com/95boi.jpg) I am having trouble finding expected values. As the question states, the expected values can be generated from null hypothesis. But how? Another confusion I'm having is identifying the procedure that I need to follow, since it's not stated in the question. Since I need to have expected values, does it mean that I need to follow goodness of fit to a Poisson distribution in order to identify independence between two variables? Thank you!


title: How to obtain a p-value for the eta measure of association?

text: How can a P-Value be obtained for an analysis of Eta association? Is there a reference for this?


title: Mean and covariance of a mixture of two distributions

text: This is a homework problem. If I could have some hints on how to solve it that would be great. Mean and Covariance of a Mixture of Two Distributions: $$ I = \left\\{ \begin{array}{l l} 1 & \quad \text{if subject in group 1}\\\ 0 & \quad \text{if subject in group 2} \end{array} \right.$$ Suppose $\mathbf{Z}=(Z_1,\ldots,Z_p)^T$ is a vector of measurements on a subject with mean and covariances: $$ E[\mathbf{Z}|I=1]=\mathbf{\mu}_1, \quad E[\mathbf{Z}|I=0]=\mathbf{\mu}_2\\\Cov[\mathbf{Z}|I=1]=\mathbf{\Sigma}_1, \quad Cov[\mathbf{Z}|I=0]=\mathbf{\Sigma}_2$$ Let $p=P(I=1)$. Find $E[\mathbf{Z}]$ and $Cov[\mathbf{Z}]$. My Attempt: (This makes sense to me) $E[\mathbf{Z}] = E[E[\mathbf{Z}|I]] = E[p\mathbf{\mu}_1+(1-p)\mathbf{\mu}_2] = p\mathbf{\mu}_1+(1-p)\mathbf{\mu}_2$. (Not sure about this one) $$ \begin{align*}Cov[\mathbf{Z}] &= Cov[E[\mathbf{Z}|I]] + E[Cov[\mathbf{Z}|I]]\\\&= Cov[p\mathbf{\mu}_1+(1-p)\mathbf{\mu}_2] + E[Cov[\mathbf{Z}|I]]\\\&= ???\end{align*}$$


title: Variational Bayes vs EP and other message-passing methods

text: I am trying to understand the difference between:   * Expectation Propagation (EP)   * Variational Bayes Wikipedia says: > Expectation Propagation differs from other Bayesian approximation approaches > such as Variational Bayesian methods. Why isn't **EP** considered a **Variational Bayes** method? Isn't EP Bayesian, and relies on message-passing to approximate the posterior? What makes a method Variational Bayes? Also, what about the following methods. Can they be considered Variational Bayes?   * Sum-product   * Mean-field methods   * Bethe-Kikuchi approximations


title: How to interpret BIC

text: I am fitting two different models to the same data. In one model, there is one free parameter for three different experimental conditions. In another model, I fit three free parameters, one for each condition. I do this for 10 subjects in a dataset. For each subject, the model with fewer free parameters has a higher BIC. But for every single subject, the difference in BIC is roughly the same (about 10). I find this very suspicious, since the BIC values themselves range from ~30 to ~1000. I have never used BIC before, and would like to say that the model with one free parameter is better.


title: Error bars for unknown distributions

text: I recently made a graph where I show the error bars for a certain number of "experiment". In another way, in my algorithm I'm minimizing the objective function so I would expect that increasing the sampling I'll get lower value of the objective function. As you can see in the graph, the second value from the left, `2.5` on the x-axis, contain only 2.5% of the configurations, so we wouldn't expect it to perform as well as if we used 100% of the configurations. I think that this is related to the **asymmetry** of the distributions. Is there any approach that can fix this problem - aka a method to compute CI for asymmetric unknown distributions? ![enter image description here](http://i.stack.imgur.com/9LGMw.png) This example should be useful to make this graph understandable! %%%%%%%%%%%%%%%%%%%% EDIT %%%%%%%%%%%%%%%%%%%%%               i = number of replicates (with different seed so different sampling every replicate)     z = objective function value     n = number of configurations     j = 1...n      Example: `n=1000`, `i=100`   1. Step 1. Analyze all the `1000` configurations and compute the minimum of `z_j`. Store it and replicate for `i`. Then compute mu and sigma of those `z_i`   2. Step 2. Analyze `50%` of the initial `1000` configuration and compute the minimum of `z_j`. Store it and replicate for i. Then compute mu and sigma of those `z_i`   3. Step 3. Analyze `10%` of the initial `1000`   4. Step 4. Analyze `5%` of the initial `1000`   5. Step 5. Analyze `2.5%` of the initial `1000` So we will have `mu_100`, `mu_50`, `mu_10`, `mu_5`, `mu_2.5` and `mu_1` and `sigma_100`, `sigma_50`,... Now I'm able to make those error bars like `mu_100 + - 2 * sigma_100`....


title: Hierarchical Bayes Normal-Normal Model

text: I have the following data for 8 runners in a 100 meter dash:               runner 1 88 91 87 82     runner 2 81 85 78 91     runner 3 75 77 83 81     runner 4 92 89 84 82     runner 5 78 79 84 92     runner 6 89 75 79 83     runner 7 91 89 92 91     runner 8 87 86 88 91      The ratings represent a performance rating and are normally distributed with unknown mean and unknown variance. Each runner can be considered as a sub- group with a mean and variance. Any guidance will be highly appreciated


title: Assessing the power of a normality test (in R)

text: I want to assess the accuracy of normality tests over different sample sizes in R (I realize that normality tests may be misleading). For example, to look at the Shapiro-Wilk test, I'm conducting the following simulation (as well as plotting the results) and would expect that as the sample size increases the probability of rejecting the null decreases:               n <- 1000     pvalue_mat <- matrix(NA, ncol = 1, nrow = n)          for(i in 10:n){         x1 <- rnorm(i, mean = 0, sd = 1)         pvalue_mat[i,] <- shapiro.test(x1)$p.value     }             plot(pvalue_mat)      My thought would be that as the sample size grows there should be a lower rejection rate, however it seems pretty uniform. I think I am misunderstanding this - any and all thoughts welcome.


title: How to compute type II error for Tukey HSD

text: I used Tukey HSD (in R) to find means that are significantly different from each other, but now I’m wondering how to estimate the type II error to have an idea how likely it is that false H0 is incorrectly accepted for pairs of objects which are not significantly different. I would like to use this to justify conclusion on non-existing differences. Unfortunately I was not able to find any helpful information on the internet (maybe because I don’t know appropriate keywords). Is there a way to compute the type II error for Tukey HSD? Or should I use a different approach to justify conclusion on non-existing differences? Hopefully someone can point me into the right direction.


title: Power of maximum likelihood parameter estimates for a linear model

text: The maximum likelihood parameter estimates for the linear model where $\Pr(Y|X\beta) \sim \mathcal{N}(0,\sigma^2)$ are: $$\hat{\beta} = (X'X)^{-1}X'Y$$ How do you compute the statistical power of the parameter estimates? It seems that there should be a closed-form, but I've searched and haven't found a proof anywhere else online.


title: Interpreting probability conditions from question

text: I've encountered this question: ![http://i.stack.imgur.com/zlE1w.jpg](http://i.stack.imgur.com/LmgEf.jpg) And got the answer here: ![enter image description here](http://i.stack.imgur.com/W2jMY.jpg) However, what I don't quite understand is how are the two conditions derived from the question in the first place.


title: Repeated measures ANOVA with unbalanced data

text: I am a master's student, and I need your help. I did experiment on one group of rats, the group was tested during 8 days, this is the first stage. The same group was also tested again but this time the test phase was 7 days. I tried to do my analysis using ANOVA repeated measures but I had a problem with the days. I want to know the difference in the 2 tests but the days are not the same number. I don't know how to handle it. If there is any way to do the analysis using repeated measures ANOVA. I tried to use repeated measures ANOVA within subjects effect. My first factor will be `test` with 2 levels (test 1 & 2), second factor is `day` with the level of? Here is the problem! If I write 8 there will be extra day in test 2 and the analysis can't run!


title: Calculation of MA model in R

text: Trying to manually calculate the fitted values from arima in R in order to understand the model. What I am finding is that even though I have not yet learned to calculate the first p values in a AR(p) or MA(p) model, even when I take values returned from arima the first few fitted values dont seem to match. What am I doing wrong? R Code:               install.packages("forecast")     library(forecast)          wineindT<-window(wineind, start=c(1987,1), end=c(1994,8))           write.table(wineindT,"clipboard",sep="\t")     mod3<-Arima(wineindT,order=c(0,0,1))     mod3$coef            write.table(mod3$coef,"clipboard",sep="\t")            write.table(forecast.Arima(mod3,1)$fitted,"clipboard",sep="\t")     forecast.Arima(mod3,1)$fitted     mod3$residuals           write.table(mod3$residuals,"clipboard",sep="\t")      Here is an Excel calculation: Given the first residual (from arima) my calculation is off for series value 2 and 3 then is correct. ![enter image description here](http://i.stack.imgur.com/bHDTo.jpg)


title: Convergence sensitive to choice of baseline category?

text: While running a simulation in R, I noticed that in fitting one particular model, R spits out a warning message, but if I simply change the baseline category in the response variable, it converges without complaining. The results produced do not differ appreciably, except that one value in the random effect correlation matrix is changed. I am curious why this occurs, and whether it suggests something interesting about the data set (the response variable is binary valued, and the explanatory variables are continuous- valued).               > library(lme4)     > fit1.MS.eE <- glmer(label ~ zSpec * zF1 + zF2 + (1 + zSpec + zF1 + zF2 | part), data = MSall.eE, family = binomial())     Warning message:     In mer_finalize(ans) : singular convergence (7)     > MSall.eE$label <- relevel(MSall.eE$label, "e")     > fit1.MS.eE <- glmer(label ~ zSpec * zF1 + zF2 + (1 + zSpec + zF1 + zF2 | part), data = MSall.eE, family = binomial())     >(no warning message)      I have uploaded the data set referred to in the code for interested persons.


title: Interpreting a weak model proceeding from there

text: I originally posted this question and since then I've posted another more specific branch-off question. I've followed the advice I received in response to my questions, and I've done some statistical analysis. Now I need help analyzing the results and determining my next move. Here's my original explanation of the situation: > My data is a list of records, each one representing an educational seminar > event. I have a continuous variable that represents the revenue brought in > by each seminar, which is the response variable in my regression. I also > have a number of categorical variables which are acting as factors/IVs. To add a little more detail, these factors include things like day of the week the seminar was held, topic, speaker, etc. My main goal is to build a model that A) can help to explain what factors most influence revenue and _how_ , and B) has some predictive power. I performed a multiple regression in R, but much to my dismay, the adjusted $R^2$ value was a mere 0.2188. I know this doesn't mean the factors have _no_ predictive power at all, but I'm wary of making any major strategic decisions based on such a poorly fit model--am I right to feel that way? I also performed an ANOVA test, and certain factors seemed to be labelled as more significant than others. My questions are:   1. How should I proceed given the weak fit of my model? Could it be that the data simply isn't helpful and that either random chance or factors I haven't considered are at play here?   2. Though ANOVA tells me which factors appear significant, it doesn't give me any idea of **how** they are significant. How can I determine what effect each significant factor has from a practical standpoint?   3. Are there any other tests that would help me understand how all these variables are related? Bear in mind that I am new to R and rusty on statistical methods (I took an intro-level course a few years ago).


title: Singularity issues in multinomial model using R

text: I am trying to develop a mode choice model (4 modes: hov, transit, bike, walk) and below are two approaches I am using. I am having problems in both   * **Approach 1** Mode choice as a function of price (cost of using the mode) `price`: generic variable   `trans1`: choice variable (0,1)   Dataset: `work`   Command:              > mode.choice <- mlogit(trans1 ~ price, work)          Error: Error in solve.default(H, g[!fixed]) :        Lapack routine dgesv: system is exactly singular      P.S: Unlike previous posts, I don’t have NA’s in my dataset.   * **Approach 2** Mode choice as a function of price and some alternative specific variables `price`: Generic Variable   `trans1`: choice variable (0,1)   `hh1`, `hh2`, `hh3`: alternative specific variables   Dataset: `work`   Command:              > mode.choice <- mlogit(trans1 ~ price | hh1 + hh2 + hh3, work)          Error: Error in solve.default(H, g[!fixed]) :        Lapack routine dgesv: system is exactly singular      I have tried different variables in both approaches but the singularity issue persists Help on any of these approaches would be greatly appreciated.


title: GEE with exchangeable working covariance  vs. assuming independence and using Huber-White standard errors?

text: I'm analyzing a dataset including 13000 students. Students are clustered into schools/grades. The ICC (intraclass correlation coefficient) shows that students in a same school are correlated. Therefore, I'd to take this clustering effect into account. One way is to run a linear regression and run the robust variance estimator on top of that to guard against getting biased estimates. **Can we take the clustering effect into account with the sandwich estimator?**. I think we cannot since the epsilon^2 matrix is still a diagonal matrix. What I think we should do is to run a GEE model with exchangeable varcov matrix for students in a same school and then we should run the robust varcov on top of the GEE model. Anyways, **assuming that sandwich estimator by itself can account for the school clustering effect** , I've fitted two separate GEE models one with exchangeable varcov matrix and the other one with the robust variance estimator (also known as Huber-White, Sandwich Estimator, or empirical variance estimator). The point is under both models I get the same estimated variance per each covariate, but my GEE exchangeable estimates leads to much larger beta estimates that are also statistically significant whereas similar beta covariates are not significant in GEE with robust varcov estimator. I wonder why it happens? I use SAS 9.3 to fit my model (proc GENMOD): Exchange:               repeated subject = SCHIID/type = EXCH;      Empirical Estimator:               repeated subject = SCHIID/ covb;


title: Independence of numerator and denominator connected to independence of obs?

text: I was wondering whether the independence of the numerator and the denominator in the ANOVA F-statistic relates to the independence of observations and how?


title: Variance of resistors in parallel

text: Suppose you have a set of resistors R, all of which are distributed with mean μ and variance σ. Consider a section of a circuit with the following layout: (r) || (r+r) || (r+r+r). The equivalent resistance of each part is r, 2r, and 3r. The variance of each section would then be $σ^2$, $2σ^2$, $3σ^2$. What is the variance in the resistance of the entire circuit? After sampling several million points, we found that the variance is approximately $.10286\sigma^2$. How would we arrive to this conclusion analytically? Edit: Resistance values are assumed to be normally distributed with some mean resistance r and variance $σ^2$.


title: Why is the lower limit of my lower 95% confidence interval negative?

text: I'm conducting some research related to the use of patient data and quality of care. I'm generating a linear regression graph between phone bill size and customer population. Assume that the intercept of this graph came up a positive value, but the 95% confidence interval for the intercept was (-50, +50).   * Since a phone bill can never be a negative value, how is it that my lower CI limit is -50?    * Does this represent an error in my calculations? (I doubt this, as I've been re-creating an old study performed by experts.)   * Or is this expected behavior?


title: Does $Y=\alpha X + \beta$ hold for multivariate gaussian density?

text: In the one-dimensional case, if $X$ is $\mathcal{N}(\mu,\sigma^2)$, then $Y =\alpha X + \beta $ is $\mathcal{N}(\alpha \mu + \beta,\alpha^2\sigma^2)$ . We can prove this using the cumulative distribution function of of $Y$ $F_Y(a) = P\\{Y \leq a\\} = P\\{\alpha X + \beta \leq a\\} = P\\{X \leq (a-\beta)/\alpha\\}$. Substituting $Y =\alpha X + \beta $ and change of variable gives us, $F_Y(a) = \int_{-\infty}^{a} \frac{1}{\sqrt{2\pi}(\alpha\sigma)} \exp \\{ \frac{-(v-(\alpha \mu + \beta))^2}{2(\alpha\sigma)^2}\\} dv $ Hence $f_Y(v) = \frac{1}{\sqrt{2\pi}(\alpha\sigma)} \exp \\{ \frac{-(v-(\alpha \mu + \beta))^2}{2(\alpha\sigma)^2}\\} $ Thus $Y$ is $\mathcal{N}(\alpha \mu + \beta, \alpha^2 \sigma^2)$. In the multivariate case, if $X$ is $\mathcal{N}(\mu,\Sigma)$ and $Y=\alpha X + \beta$, is $Y \sim \mathcal{N}(\alpha \mu + \beta,\alpha^2\Sigma)$? If so, how do we prove it?


title: How to convert multiple ranking scores into a probability distribution?

text: I would like to create a topic distribution for a document. The current model I am trying to implement is: for each sentence in the document, I am getting a topic assignment with a score, e.g. "1st sentence is about Microsoft with a relevance score of 0.4". I repeat this for each sentence, and at the end I have relevance scores with the topics like the following: > 1st sentence: microsoft, score 0.4 > > 2nd sentence: apple, score: 0.1 > > 3rd sentence: android, score: 0.5 > > ... Now, I would like to convert these scores into 1 big probability distribution that will represent the whole document. Is there a known technique to do this? If so, what is the best way to do this? Note: I know this is a very naive topic modelling, but I am currently interested in combining the scores into prob. distributions.


title: Few random variables cannot influence $n$ independent others too much?

text: I have $n$ standard normal and independent random variables $X_i$ (In reality I have a large known number of them, but let's just say I have $n$). In my experiment I want to on average get exactly 3 random variables $X_i$ under a threshold $c$. To get that, I can compute $c$ having that property easily, because the average number of $X_i$ that are under a threshold $c$ is $n \Phi(c)$ where $\Phi$ is the cdf of the standard normal distribution. So I choose $c = \Phi^{-1}(3/n)$ in this case. (Which is a negative number for large $n$.) But unfortunately I already know the value of two other standard normal random variables $Y$ and $Z$ which may depend on each other and on any number of the other $X_i$. So my question is: if I know that $Y$ and $Z$ are under the threshold $c=\Phi^{-1}(3/n)$, is it then still true that on average at most a constant number of the other random variables $X_i$ are under the threshold $c$? So by knowing that $Y$ and $Z$ are under the threshold, they can't suddenly make many of the other random variables go under it too. I am almost certain that they can't, but I don't know how to prove it. Any hints are welcome. Or books where you think this might be in.


title: Estimating right-censored data

text: I am VERY new to stats. I have a large amount of life-time data (delay in arrival since start of experiment) from repeat experiments. Some data is missing, but essentially represents a delay longer than the experiment time (right-censored data). I want to estimate this data. I could use the experiment time for each missing data, but I understand this is not useful above a certain data size. I understand there are statistical strategies to estimate using the distribution of the non-missing values, the fact that the missing data is right-censored and some kind of maximum likelihood equation to find the parameters (missing data) that would have the highest probability of giving the non-missing data, given the distribution. However, this is not the larger part of my work, and I am really looking to use some kind of tool, rather than endless manual crunching. Or... If there were available a very clear pseudo code algorithm, I could easily port it to C/C++, because the original data is generated in my C++ application. The data is later offloaded into Excel and so alternatively achieving the same thing in Excel that would be nearly as useful, but not quite. Any suggestions would be appreciated.


title: Centroid matching problem

text: For a Dataset $D$, we have gold standard centroids say $c_1, c_2, \cdots, c_n$. Now if we run k-means algorithm on $D$ with input $n$, we get k-means centroid $k_1, k_2, \cdots, k_n$. I just wanted to know, is there any algorithm/heuristic to match the centroids between $k_i$ and $c_j$ where $i, j= 1, \cdots, n$ (One to one mapping between $k$'s and $c$'s) I tried to calculate the pairwise distance between $k_p$ and $c_j,\; j= 1, \cdots, n$, and match $k_p$ to $c_r$ where the distance between them is minimum. But in this case $k_p$ and $k_q$ are assigned to $c_r$, which we dont need.


title: Is it statistically correct and safe to calculate Kappa 2 test means of combinations of 2 against the value of a Fleiss' Kappa?

text: For my problem I have a group of `p participants` and `s subjects`. I was suggested to calculate the Kappa Cohen Test on combinations of two for `p` and `s`, then calculate the `mean` and compare against the uniquely calculated value of Fleiss` kappa for the same set. My question is, whetheror not it is correct/safe to calculate a mean over the Cohen's kappa test? More specifically, I am using the Kappa2 and Fleiss' kappa with the flag set to FALSE from the irr R package. Setting the flag to false means that a variant algorithm of Kappa Fleiss is used so that it can be reduced to Kappa Cohen when the population has p=2 s=2.


title: Which test do i choose to conduct on my Likert scale data?

text: I have been doing plenty of research into which test I should use on my Likert scale data and i am becoming more confused as there are so many people saying different things. I am currently doing my dissertation on To what extent does music affect aggression whilst playing violent and nonviolent video games? I have use three different methods to collect data from 30 participants:   1. An aggression scale before and after the participants played both violent and non-violent video games to measure the changes in aggression whilst listening to different genres of music   2. Story completion tasks which give the participant a scenario and asks for ten responses i.e i would be angry, I'd hit the driver etc after playing both violent and non-violent video games   3. A questionnaire to see what peoples' opinions are on violence in video games using a likert scale ( 5 point scale strongly agree - strongly disagree) and also open ended questions on preferences to music etc. I have already completed my analysis for the aggression scales (1) but i am completely stuck on which test i should choose when trying to gather data from my likert scales. the questions were such as:   * Video games are becoming too violent    * Violence in video games are influencing people to become more aggressive    * Video games are unsociable    * Video games teach and encourage people's aggressive behaviour   * Video games help people to forget about everyday life worries   * Video games encourage people to be more aggressive in real life situations which they then score on a scale of Strongly Agree, Agree, Neither Agree nor Disagree, Disagree, and Strongly Disagree for all statements (13 all together). I want to know whether people believe games cause real life violence or whether they think games are too violent. Thanks in advance :)


title: Cross variogram with a moving window

text: I need to generate cross variograms of images using moving windows. For that I use the following equation: $$ \gamma_{jk}(h)=\frac{1}{2n(h)}\sum_{i=1}^{n(h)}\Big\\{\big[dn_j(x_i)-dn_j(x_i+h)\big]\cdot\big[dn_k(x_i)-dn_k(x_i+h)\big]\Big\\} $$ The first part stands for one band(j) and next part of band k. To illustrate with sample matrices,               j =  1     2     3     4          5     6     7     8          9    10    11    12         13    14    15    16          k = 17    18    19    20         21    22    23    24         25    26    27    28         29    30    31    32      In actual case I am using 7 X 7 windows for large satellite images. I also had to generate variograms from images for this work. For generation of variograms I had to consider only one band of data. For that case I used `nlfilter` for moving window and created a function to select and calculate values. But for cross variograms, I am not able to decide upon what function to use. For this case the calculations will go like this: > (1 - 2)(17 - 18) + (2 - 3)(18 - 19) and so on.


title: Percentile calculation based on a known percentile table

text: I am a very newbie of R. Please help me out of this simple question. Thanks. Now, I have two data about people's weights. One is from normal patients, and another is about patients with a specific disease. Let's say:               weight_Normals <- c(170, 140, 122, 143, 156, 122, 189, 200, 201, 133, 155, .....)     weight_Disease <- c(220, 145, 167, 177, 190, 130, 169, 100, 230, 210, 222, .....)      Now, I can compute the percentile from 0.01 to 0.99 of weight_Normals. However, how should I use this percentile results from normal patient to map into the patients with diasease? In other words, I would like to calculate the percentile of these patients with disease based on the percentiles of normal patients. Please let me know if my question is not clear enough. Appreciate your time and help.