title: Conditional distribution of uniform random variable distributed over (0,1)

text: Let $U$ be a random variable uniformly distributed over (0,1). Compute the conditional distribution of U given that $U>a$ The solution says: $P(U > s | U > t) = \frac{P(U > s)}{P(U > t)}$ Here I'm confused why the conditional probability translates into the fraction on the right hand side. $P(U > s | U > t) = \frac{P(U > s)}{P(U > t)}$ $=\frac{1-s}{1-t}$ where $t \lt s \lt 1$ The distribution function of a uniform variable $P(U \leq x) = \frac{x-a}{b-a}$, so in this case, $P(U > s)=\frac{s-t}{1-t}$ because $t \lt s \lt 1$. $P(U > t)=\frac{t-0}{1}=t$ so $\frac{P(U > s)}{P(U > t)}=\frac{\frac{s-t}{1-t}}{t}$ however this is different from the answer in the solutions, which is $\frac{1-s}{1-t}$. Can someone please explain why my answer is incorrect?


title: Survival analyses for complicated models

text: I've just began looking into survival analyses, and I'm having some difficulty interpreting the results. I have a model linear regression model where time ~ x1 * x2 * x3... x9, this model seems to fit well (adjusted R${^2}$ of ~0.9) however it was suggested that survival analyses may be a better direction. I have tried using the R tool for this, but I can't interpret the data. I do have some trails where they timed out (after 1 hour) so I set up my censor column for that. I then tried plotting the full data set (~200 trials), and got a nice downward stepped curve. ![Full data set](http://i.stack.imgur.com/JVEiV.png), however this doesnt seem to tell me anything, the model used was just survfit(my.surv ~ 1). If I replace the ~ 1 with one of the variables, for example survfit(mysurv ~ x1), I get a total mess. Lines all over the place (mostly vertical). So, my questions are these: 1) Is it possible to gain useful information, other than saying that 80% of the trials finished before 2000 seconds, which is what I think this plot shows 2) Is it possible to include factors in the model, for example time ~ x1 + x2, etc. If not how does one ascertain which variables are important by this method. 3) How can I use this to make predictions about new data, not included in this model?


title: Segmented regression in R: interpreting the other coefficients

text: I am running segmented regression using the R package 'segmented'. The original binomial logistic regression has two coefficients, approach_km (continuous), and sea (dichotomous) that explain the probability of the dependent variable occurring. I am interested in investigating if there is a break-point in the regression slope for approach_km. Here are the coefficient estimates from the original non-segmented glm               Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) -1.15563    0.17352  -6.660 2.74e-11 ***     approach_km -0.03915    0.01140  -3.432 0.000598 ***     sea2         0.34501    0.20515   1.682 0.092615 .       Here are the estimates after running the model through the segmented function               Estimated Break-Point(s):        Est. St.Err       3.794  1.134           t value for the gap-variable(s) V:  0           Meaningful coefficients of the linear terms:                    Estimate Std. Error z value Pr(>|z|)       (Intercept)    -0.08233    0.32692  -0.252    0.801       approach_km    -0.49096    0.19776  -2.483    0.013 *     sea2            0.27182    0.21141   1.286    0.199       U1.approach_km  0.49085    0.19817   2.477       NA        A Davies' test indicates that the differences in the two slope estimates is significant. The question I have is how do I interpret the coefficient sea? I see that the estimate (and the p value) changes in the segmented output. Does this mean that this is the best estimate for the value of the coefficient for both of the estimated slopes? Since the sea coefficient is not significant we can assume that this variable will not statistically alter the overall effect or interpretation of the model but I am interested in how I would interpret this result if it was in fact significant. Thanks.


title: Spline – basis functions

text: I am trying to wrap my head around splines and the concept of basis functions using the _Elements of Statistical Learning_. I understand that the goal is to find polynomials that are continuous at first and second derivatives. However, following the picture below, I don't understand whether   * a) the spline consists of a different cubic function $(a+bx+cx^2+dx^3)$ in each of the three regions, or    * b) whether the spline is the linear addition of the 6 basis functions per below across the entire domain, or    * c) whether there are 6 basis functions with different parameters in each of the 3 regions (hence 18 different functions). Much appreciated... ![enter image description here](http://i.stack.imgur.com/6Audn.png)


title: Accurately generating variates from discrete power law distribution

text: What are the best methods to _accurately_ generate random integers distributed according to a power law? The probability of getting $k$ ($k=1,2,\ldots$) should be equal to $p_k = k^{-\gamma} / \zeta(\gamma)$ and the method should work well for any $\gamma > 1$. I can see two naive approaches:   1. Compute $p_k$ up to some large $k_\text{max}$ so that $\sum_{k=1}^{k_\text{max}}$ is "close enough" to 1, then generate integers according to these probabilities. This just won't work if $\gamma$ is close to 1 as $k_\text{max}$ would need to be huge.   2. Draw real numbers from a _continuous_ power law distribution (an easier problem that I know how to solve) and round them to integers in some way. It is possible to analytically compute the precise probability of obtaining each integer with the above method. I could use rejection to correct these to $p_k$ (which can also be computed provided I can evaluate the $\zeta$ function). (This would be a bit hairy as I'd have to round in a way that I obtain integers with higher probability than $p_k$ for $k$ greater than some small value, and handle $k$ less than that separately.) Is there a better method that is also accurate (not approximate)?


title: Probability - overlapping normal curves

text: Two cyclists are racing over 50 (for example) km. Cyclist A has an estimated speed of X km per hour, cyclist B has an estimated speed of Y km per hour. Both of these estimated speeds are the means of two normal curves of the true actual speed. The two normal curves overlap. If I know the estimated speeds of both A and B, and the respective standard deviations of their normal curves, how do I calculate the probability of A winning and of B winning? What about when there are n cyclists, not just two? And, if I only know the probabilities of A winning and B winning, and I guess the standard deviations of their normal curves, how can I calculate their relative speeds? And for n cyclists? Thanks


title: How can I use R to build an ordered logistic regression model and to measures its performance?

text: I want to use R to build an ordered logistic regression model and to measures its performance. My goal is to predict whether an athlete will get gold, silver or bronze. I started with logistic regression where the number of available categories is two, e.g. won and lost.               str(DATA$dependent_variable)     Factor w/ 2 levels "won","lost",...      So far, I'm satisfied with my code for logistic regression:               require(caret)     DATA <- read.csv("data.csv", header = T, row.names = NULL)     set.seed(4444)     trainIndicator = rbinom(nrow(DATA),size=1,prob=0.5)     DATA.train = DATA[trainIndicator==1,]     DATA.test = DATA[trainIndicator==0,]          model <- train(dependent_variable ~ predictor1 + predictor2 + predictor3, data = DATA.train, method = "glmStepAIC", direction = "both", family = binomial(link = logit), trace = F)          DATA.predicted <- predict(model, DATA.test)     print(confusionMatrix(DATA.predicted, DATA.test$dependent_variable))      The last step shows me performance characteristics such as accuracy, sensitivity, specificity, etc. How do I build an ordered logistic regression model in R and how do I measure the performance of an ordered logistic regression model when the dependent variable has three categories, e.g. gold, silver and bronze?               str(DATA$dependent_variable)     Factor w/ 3 levels "gold","silver", "bronze",...


title: Steps done in factor analysis compared to steps done in PCA

text: I know basically how to express PCA (Principal component analysis) mathematicaly, but I would like to know steps that should be used for factor analysis.               Let us consider some matrix, for instance              A=[3 1 -1;2 4 0;4 -2 -5;11 22 20]              A =                   3     1    -1              2     4     0              4    -2    -5             11    22    20          i have calculated correlation matrix              B=corr(A)              B =                  1.0000    0.9087    0.9250             0.9087    1.0000    0.9970             0.9250    0.9970    1.0000          now i have done eigenvalue decomposition              [V ,D]=eig(B)              V =                  0.5662    0.8209   -0.0740             0.5812   -0.4613   -0.6703             0.5844   -0.3366    0.7383                   D =                  2.8877         0         0                  0    0.1101         0                  0         0    0.0022      General idea behind the PCA is to choose significant components, form new matrix which have columns eigenvectors, then we need to project original matrix (in PCA it is zero centered). But in factor analysis, for instance we should choose components that have higher that $1$ singular value, also then we are using rotation of factors, please tell me how it is done? For instance in this case, so please help me **understand factor analysis steps relative PCA steps**.


title: extractSeqFromW function

text: I'm using R and the library TraMineR to do sequence analysis on social science data. I need to use the extractSeqFromW function, but I have problems opening the file even if they are exactly in that directory.               source("extractSeqFromW.R")          sla <- extractSeqFromW("Users/emanuelastruffolino/Desktop/SequenceCourse/Assignment6","Users/emanuelastruffolino/Desktop/SequenceCourse/Assignment6",  pvarseq="P$$A06",            use.value.labels=FALSE, maxMissing=5)      This is the error message               Error in read.spss(paste(datadir, "MASTER.sav", sep = ""), to.data.frame = TRUE) :        unable to open file: 'No such file or directory'      What's wrong? Thank you!


title: Using glmer to replicate result from lmer for mulitlevel modelling in R

text: Having read from the internet that `glmer` and `lmer` are similar thing (correct me if I am wrong), just `glmer` need to state the family, I wonder if I state the family for `glmer` as "gaussian", can I replicate the same result as `lmer` does? > mlm.3 <- > lmer(formula=dv~floor*collection_date+(floor*collection_date|area_desc),data=master.df,REML=FALSE) > mlm.3.mod > <-glmer(formula=dv~floor*collection_date+(floor*collection_date|area_desc),data=master.df,family="gaussian") I tried but was unable to get it to work, so I wonder if I have some misconception here.


title: Not sure how to work out SSG and SSE in ANOVA table

text: > **Possible Duplicate:**   >  Not sure how to work out SSG and SSE on ANOVA table I have a table that is: $$ \begin{matrix} & \mathrm{Tyre\, A} & \mathrm{Tyre\, B} & \mathrm{Tyre\, C} & \mathrm{Total} \\ \mathrm{number\, of\, tyres\, tested} & 25 & 32 & 40 & 97 \\ \mathrm{mean\, response*} & 18.46 & 19.56 & 16.47 & 18.00 \\ \mathrm{sample\, variance**} & 14.42 & 13.66 & 26.77 & 20.73 \end{matrix} $$ *mean was calculated as $bar{y} = \frac{1}{n} \sum_{i=1}^{n}y_i$. **variance was calculated as $ \frac{1}{n-1} \sum_{i = 1}^{n} (y_i - \bar{y})^2 $ Using this information, I need to create an ANOVA table. The degrees of freedom were easy (I got 2, 94 and 96). But I am now stuck on the SS columns of the table. For SST I did: $(N - 1)S_T^2 $ and used 20.73 as the total variance and this gave me SST = 1990.08. Is that correct? Also, I am stuck on how to work out SSG or SSE. I know if I can work one out, then I'll have the other but how do I do this?


title: Finding anomalies in industrial part failures

text: I have a data set that is composed of the age an industrial machine part is when it fails (x variable) and the part's width in nanometers (y variable) at the time when it fails. The data looks something like this:               age (days)     width (nm)     10             100     18             88     36             45     ...            ...      I have millions of data points similar to the data above and graphically I see a nice non-linear, monotonic correlation between x and y. For a new data sample, I'd like to derive some measure as to whether this measurement was 'likely' to come from my population. The measure could be a goodness of fit statistic or probability perhaps. My question is, what is the best way to get a metric about whether a new sample was likely to come from my population? For example, I'd like to know what the likelihood or probability of drawing a (x,y) sample such as (1000, 100) from my distribution. Based on the plot, this sample is an extreme outlier, but I want a metric that quantifies this. Similarity, (10, 100) is a very frequent sample in my population and is very likely from my population, I'd like to be able to measure the same thing for this sample as well. Thanks in advance for the help.


title: Periods in history of statistics

text: The history of many fields of science can be divided into a small number of time intervals that often begin with some important discovery. But I have never seen something similar in timeline of statistics. Obviously, there are some important dates that can be considered as starting points of a new period (Pascal+Fermat, Bayes, Pearson, Tukey,..). Can we at least very roughly divide history of statistics into small number of periods? Note that the only similar question to this is related to only famous statisticians, not to periods in history.


title: In case of multi-time series data, does the significant linear model gives there is something more?

text: I have total of 9 time series data. I'm basically using 8 variables as my X and 1 variable as Y. I've tried VAR model today but I don't see ANYTHING. However, when I do linear regression, all of the variables are significant, and p value of the model is less than 0.01 %. Is this result indicates something?


title: Model comparision

text: I am comparing supervised vs unsupervised method. I have computed entropy from the confusion matrix of the 10-fold crossvalidated output. I have also applied external criteria to unsupervised methods and computed entropy from it. Next checked its significance. I was wondering when I compare supervised entropy vs unsupervised entropy should I do a significance testing for the difference in entropy between them or just describe the results as the method with lowest entropy is the best method. I will highly appreciate any help.


title: Identifying volatile "moves" in a price timeseries

text: Necessary advance apologies: I am not a statistician, so I hope I am able to make my question as clear as possible. I have a price timeseries... say, something like this, but much more granular: ![enter image description here](http://i.stack.imgur.com/RNJ3E.jpg) And I'm looking for a simple way to algorithmically identify volatile "moves", such as the large spike around "04 Aug" in the above graph. This is possible a repost of Identify different periods of variance in a time series , but I'm less interested in "periods" of volatility, more interested in the actual "moves" in one direction. So I want to identify peak-to-trough/trough-to-peak, but only if it's volatile, if that makes sense.


title: Multivariate regression

text: I have a multivariate problem. I have energy usage, max. temprature, and days. I want to predict load which is energy usage for various consumers. LM does not work as it is not linear. I am planning to use multivariate ARIMA. I want to implement it in R. I will appreciate if you could suggest some good pointers.


title: Good libraries for working with probabilistic graphical models?

text: Could someone recommend some well-maintained and up-to-date libraries for working with probabilistic graphical models? I noticed that there are some libraries for R listed here and one for C++, but are there some other good libraries in C++ or Python?


title: When do Poisson and negative binomial regressions fit the same coefficients?

text: I’ve noticed that in R, Poisson and negative binomial (NB) regressions always seem to fit the same coefficients for categorical, but not continuous, predictors. For example, here's a regression with a categorical predictor:               data(warpbreaks)     library(MASS)          rs1 = glm(breaks ~ tension, data=warpbreaks, family="poisson")     rs2 = glm.nb(breaks ~ tension, data=warpbreaks)          #compare coefficients     cbind("Poisson"=coef(rs1), "NB"=coef(rs2))      ![enter image description here](http://i.stack.imgur.com/cJkk6.png) Here is an example with a continuous predictor, where the Poisson and NB fit different coefficients:               data(cars)     rs1 = glm(dist ~ speed, data=cars, family="poisson")     rs2 = glm.nb(dist ~ speed, data=cars)          #compare coefficients     cbind("Poisson"=coef(rs1), "NB"=coef(rs2))      ![enter image description here](http://i.stack.imgur.com/G2Q9i.png) (Of course these aren't count data, and the models aren't meaningful...) Then I recode the predictor into a factor, and the two models again fit the same coefficients:               library(Hmisc)     speedCat = cut2(cars$speed, g=5)      #you can change g to get a different number of bins          rs1 = glm(cars$dist ~ speedCat, family="poisson")     rs2 = glm.nb(cars$dist ~ speedCat)          #compare coefficients     cbind("Poisson"=coef(rs1), "NB"=coef(rs2))      ![enter image description here](http://i.stack.imgur.com/25sbs.png) However, Joseph Hilbe’s _Negative Binomial Regression_ gives an example (6.3, pg 118-119) where a categorical predictor, sex, is fit with slightly different coefficients by the Poisson ($b=0.883$) and NB ($b=0.881$). He says: “The incidence rate ratios between the Poisson and NB models are quite similar. This is not surprising given the proximity of $\alpha$ [corresponding to $1/\theta$ in R] to zero.” I understand this, but in the above examples, `summary(rs2)` tells us that $\theta$ is estimated at 9.16 and 7.93 respectively. So why are the coefficients exactly the same? And why only for categorical predictors? * * * ## Edit #1 Here is an example with two non-orthogonal predictors. Indeed, the coefficients are no longer the same:               data(cars)          #make random categorical predictor     set.seed(1); randomCats1 = sample( c("A","B","C"), length(cars$dist), replace=T)     set.seed(2); randomCats2 = sample( c("C","D","E"), length(cars$dist), replace=T)          rs1 = glm(dist ~ randomCats1 + randomCats2, data=cars, family="poisson")     rs2 = glm.nb(dist ~ randomCats1 + randomCats2, data=cars)          #compare coefficients     cbind("Poisson"=coef(rs1), "NB"=coef(rs2))      ![enter image description here](http://i.stack.imgur.com/y7S7K.png) And, including another predictor causes the models to fit different coefficients even when the new predictor is continuous. So, it is something to do with the orthogonality of the dummy variables I created in my original example?               rs1 = glm(dist ~ randomCats1 + speed, data=cars, family="poisson")     rs2 = glm.nb(dist ~ randomCats1 + speed, data=cars)          #compare coefficients     cbind("Poisson"=coef(rs1), "NB"=coef(rs2))      ![enter image description here](http://i.stack.imgur.com/TxyuH.png)


title: How to handle gaps in a time series when doing GAMM?

text: I want to apply a GAMM with R to this time series but I am not sure how to handle the station P18, as shown in the figure below. ![enter image description here](http://i.stack.imgur.com/4uSak.png) If I shrink the dataset to the point where P18 ends (i.e. left side of the dashed line) I would loose many reversing trends in P07,P14,P22,P2O and so on. The question is thus: can I use the whole dataset including P18? I am not sure what the GAMM is doing with such a gap. what would be if there were more gaps? any idea is welcome!


title: Comparing fatality rates

text: I have case fatality rates (deaths per 100 cases) for 2 different states receiving different treatments for 17 years. What is the best statistical method to compare them? Relative risk, odds ratio, plain time series analysis, ...? The data is like this:               Year   St.1 Cases  St.1 Deaths St.1 CFR St.2 Cases St.2 Deaths  State2 CFR       1994    1836        383         20.86     583      121          20.75     1995    1246        257         20.63     1126     227          20.16     1996    1450        263         18.14     896      179          19.98     1997    2953        407         13.78     351      76           21.65     1998    1161        149         12.83     1061     195          18.55     1999    2924        434         14.84     1371     275          20.06     2000    1729        169         9.77      1170     253          21.62     2001    1888        275         14.57     1005     199          19.80       2002    919         178         19.37     604      133          22.02        2003    865         142         16.42     1124     237          21.09     2004    1543        131         8.49      1030     228          22.14     2005    2887        336         11.64     6061     1500         24.75     2006    1484        108         7.28      2320     528          22.76     2007    1592        75          4.71      3024     645          21.33     2008    1920        53          2.76      3012     537          17.83     2009    1477        40          2.71      3073     556          18.09     2010    1534        26          1.69      3540     494          13.95


title: ar() time series function in R, manually checking the residuals/predicted values

text: I am using the ar() function to fit an AR model to some data, and this object will return the in sample residuals. I also know the syntax for how to get the corresponding predicted values, but I want to compute these predicted values manually (just to check my own understanding) for a simple AR(1) example. My problem is that my manually computed residuals (based on my manually computed predictions) do not match the in sample residuals stored in the ar object (well the 1st residual does match, but not the rest). From the documentation, I see that               x[t]  = m + a[1]*(x[t-1] - m)       where m is the sample mean of the series. Here is an example of what I am doing manually.               # Create some true AR(1) data     set.seed(123)     x = w = rnorm(30) + 3 ; for (t in 2:30) x[t] = .60*x[t-1] + w[t]     # Fit an ar model     x.model = ar(x) # coefficient is .49, mean value of x is 6.98     # Manually create predictions     x.MyPred = rep(0,30) ; x.MyPred[1] = x[1]     for (t in 2:30) x.MyPred[t] = 6.984234 + .4988327*(x.MyPred[t-1] - 6.984234)     MyResid = x - x.MyPred     cbind(MyResid, x.model$res) # does not match      And interestingly, the first residual (at observation 2) does match, but the rest do not. Thanks in advance. **Update** The answer given below is basically highlighting the difference between so called static forecasting and dynamic forecasting, here are a few more details. The two possible choices to make the fitted values are               # Method 1:     for (t in 2:30) x.MyPred[t] = 6.984234 + .4988327*(x.MyPred[t-1] - 6.984234)          # Method 2:     for (t in 2:30) x.MyPred[t] = 6.984234 + .4988327*(x[t-1] - 6.984234)      Method 1 is taking the forecasted value from the previous step and using this to plug in for the lagged values on the RHS. This could feasibly be used to create out of sample predictions forever. Also, depending on what time point you start to make predictions, the resulting predictions for a given time period can be different. Method 2 is taking the actual known historical values to plug in for the lagged values on the RHS. It will never be able to forecast more than 1 step out of sample. Also, the forecast values will always be the same no matter where you start. Both of these will give the same 1 step ahead forecast value. Method 2 is what produces the residuals returned in the ar() object. The Eviews software documentation has a good discussion of this. This R documentation was problematic in not specifying which of the x values are fitted values versus known historical values. A better expression to show in the documentation would be               \hat{x}[t]  = m + a[1]*(x[t-1] - m)      and, dare I say, possibly a few sentences on this very topic. But good documentation is hardly something to expect in R.


title: General questions about a mixed repeated measures model

text: I am using a mixed repeated measures with an ARMA(1,1) variance covariance structure. I have a basic understanding of why I am using this model etc but I need some help on articulating a few things:   1. How does this model correct for the variance in time between each time point? Or does it ignore the problem? In which case then why is it acceptable here?    2. What happens to variables that are correlated? When it is assigning variance to one variable or the other? How does this analysis do this?   3. Why this over repeated measure? Less error? More power? How much is imputed?  I have formulated some thought/answers, but I am curious to hear from others who might actually have used this model more than me. Thank so much


title: Show that $E(XY) = σ_{XY} + μ_Xμ_Y$?

text: Let X and Y be two random variables with means $μ_X$ and $μ_Y$ , variances $σ^2_Y$, $σ^2_X$ and covariance $σ_{XY}$ . Show that $E(XY) = σ_{XY} + μ_Xμ_Y$? I tried the Covariance formula with: $$Cov(X,Y) = E(X − E(X))(Y − E(Y ))$$ $$\int_{-\infty}^{\infty}{(x − E(X))(y − E(Y ))f(x, y)dxdy}=$$ $$=E(XY ) − E(X)E(Y )$$ $$σ_{XY}=E(XY) - μ_Xμ_Y$$ $$E(XY) = σ_{XY} + μ_Xμ_Y$$ Is this assumption right, because I do not need the two variances in my approach?


title: Plotting averages when there are missing values

text: I have count data for a number of subjects in different groups I would like to compare. The averages of the cumulative sums are shown in the figure below. As you can see the Red group has a spike at the end. This is because one subject in that group had fewer trials. This subject had the lowest cumulative number of counts for the previous trials and was pulling down the average. ![enter image description here](http://i.stack.imgur.com/BSslX.png) How would you plot and interpret data like this? I wish to avoid imputation or other techniques "fancier" than calculating simple statistics, ideally I would plot only the data itself. I already have a plot containing each individual, but also want a group comparison. **Edit:** The y axis is cumulative # of successes while learning a new task at each block of trials. The x axis should probably read "Block #" rather than trial number. The successes were followed by a reward, so cumulative number of successes is also a measure of total reward received over the course of the study, which may be the outcome the subjects were attempting to attain. Alternatively the learning can be measured as the number of successes per block over time. Some subjects may learn slower but end up with higher final performance and thus appear "better", however subjects that used a "quick and dirty" strategy may be maximizing the total reward attained despite their apparent lower skill level.


title: Finding a proper t-distribution given a constraint about the weight of its tail

text: I wish to find the parameter $\nu>2$ for the student t-distribution such that the following constraint will hold: if $F_T$ is the CDF of that distribution and $F_N$ is the CDF of the normal distribution, then $$\frac{F_T\left(-2\sqrt{\frac{\nu}{\nu-2}}\right)}{F_N(-2)}=c$$ for some pre- determined $c$. The motivation is as follows: the ratio between the CDFs is some measurement of weight of tails. Since both distributions are symmetric, it should suffice to check the left tail. The choice of $-2$ is arbitrary, and due to previous analysis I've made. The reason $\sqrt{\frac{\nu}{\nu-2}}$ appears there is since this is the standard deviation of the t-distribution with parameter $\nu$, and the "definition" of a left tail here is whatever appears beyond 2 standard deviations. I assume one may not hope to find $\nu$ for any given $c$. For $c=1$ the answer must be $\nu=\infty$; I believe it is impossible for $c<1$. I don't know about any upper bound, but I believe there is such. Currently I am interested in $c$ values which are in the interval $(1,2)$.


title: How to test the variance in timeseries?

text: I have a doubt regarding the variance, I try to explain It with an example. I have two vectors, like:               a <- c(1:10)     1  2  3  4  5  6  7  8  9 10          b <- c(10:1)     10  9  8  7  6  5  4  3  2  1      the variance is obviouly the same:               > var(a)     9.166667     > var(b)     9.166667      Ok, I need to test if the variances are similar, and for this test I use **var.test()**. The problem is that the variances are equals, OK! but the follow a different direction, the first move from 1 to 10 and the second from 10 to 1. SO the variances are the same and the test pass successfully(obviously), but I need also check the direction, so:   1. Are the variances similar? Ok...   2. Are the variances _(I know 'variances' here is wrong but try to understand what I mean reading the example above)_ moving in the same direction?  With the same direction I mean, the variance is equal(similar) BUT are they UP/DOWN togheter? I need to do those checks because I'm analyzing two financial lists of prices, and I need to know if the variance between their returns is constant and on the same direction. How Can I do? Thanks!


title: Optimum linear combination of predictors that maximizes a quantity

text: I use 48 energy functions to score protein-ligand interactions. I have a dataset of protein-ligands for which I know the experimental binding energy, so I can compare it with the score assigned from each energy function. With those data I can do step-wise multiple regression to find a linear combination of the following form: Pred_Energy = a*EFA + b*EFB + c*EFC +d where, EFA,EFB,EFC are 3 scoring functions from the 48, whose linear combination given above produces energies (Pred_Energy) that correlate the most to the experimental binding energies (Exp_Energy). Now from Exp_Energy I can rank the dataset of ligands according to their binding affinity, and I can compare their ranking with the ranking obtained from Pred_Energy by estimating a statistic like Kendall's W. So what I seek is a method that produces that linear combination of energy functions that maximizes Kendall's W. I would be grateful if you could give me some hints about how to do that, algorithms, keywords, anything...


title: Showing that a certain random variable has variance $1$

text: I am calculating this example at the moment. Are my calculations right? The random variable $Y$ has a mean of 15 and a variance of 9. $Z = 1/3 (Y − 15)$. Show that $Z$ has variance $1$. $Z = \frac13 (Y − 15)$ $, \newcommand{\stdev}{\mathrm{stdev}}\stdev(Y)=3$ $$\stdev(Z)=\frac13 (Y − 15) = \frac13 (Y)$$ $$\stdev(Y)=3$$ $$\stdev(Z)=\frac13 Y = \frac13 3 = 1$$ $$\mathrm{var}(Z) = 1^2 = 1$$ Is my argumentation right?


title: Can log-likelihood function calculated value (M-step) be smaller after 1 EM-iteration?

text: I am applying a MAP log-likelihood approach in order to fit a Markov mixture model, where objective function to be maximized is given by the formula: $$ L(X|\Theta _K)=\sum_{i=1}^{n}f(X_i|\Theta_K)+\sum_{j=1}^{K}\sum_{n=0}^{M}\log p(\theta_n^{j}|a_n^{j}) $$ where the second argument is a sum of Dirichlet priors ( I am using the formula given by Wikipedia) and the first argument is the sum of log- likelihood across all sequences and all components. At this point I have achieved a lot in implementing the algorithm in `R`, thanks to answers to my previously posted questions related to topic. At this stage my question is - after performing one step of expectation- maximization algorith ( with 2 components to start with), my value of $L(X|\Theta_K)$ became much smaller than it was ( from -2200 to -8000). I believe that my code is correct and do not understand why this could be happening ( the next 2 steps show steady increase). Can there be fluctuation of the algorithm in the beginning? There are 2 possible issues, however I cannot pinpoint if they are indeed the causing this: underflow problem(some of the values in my transition matrix and the resulting likelihood values, can be so negligibly small that in certain calculations `R`, in which I am working, rounds them up and thus I end up with incorrect calculations; or meaningless Dirichlet priors ( e.g. sometimes the priors are larger than 1). Additionally: the resulting from M-step multinomial distributions for each row of the transition matrix and start probabilities vector sum to 1 for each of the components. The posterior conditional probabilities of the hidden variables ( components ) also seem to make sense. The example of conditional posteriors of hidden variables for first ([[1]]) and second [[2]] component for each of 50 sequences are below ( results of first iteration E-step):                   [[1]]          [1] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00   1.520433e-10 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00         [11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00  1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00         [21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00         [31] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 2.726330e-02 1.000000e+00 1.000000e+00         [41] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.979822e-01 1.000000e+00 1.000000e+00               [[2]]          [1]  7.324479e-65  2.462187e-97  4.146568e-32  3.498317e-97 2.135013e-274  1.000000e+00  7.731884e-47 2.068553e-264 8.497501e-260          [10] 4.356689e-271 2.983088e-213 7.485688e-110 4.556750e-287 1.360219e- 173  2.340220e-45  2.609916e-59  6.057344e-59 1.286382e-185          [19]  3.879706e-80 8.488843e-188  1.881308e-14 3.098226e-290 1.681928e-290 1.168018e-211 2.123491e-292 5.767748e-177 9.232827e-198         [28] 1.120970e-159 1.397181e-257  4.078388e-48 1.524531e-247  1.000000e+00 1.833904e-252 1.452165e-263 1.878481e-111 3.379251e-178         [37] 1.823290e-247  9.727367e-01 8.553065e-238  2.748773e-32 4.602824e-138  1.212533e-93 1.744806e-271 2.677587e-292 7.822883e-131          [46] 2.504779e-111  1.775144e-16  8.020178e-01  1.301381e-63 4.437099e-114


title: Case-control study and Logistic regression

text: Suppose we have case-control data, where cases have some disease ($Y$) and controls don't and we are interested in the association of some other variable(s) ($X$). I know that in this scenario we cannot use the disease as the response variable because of the experimental design (the marginal distribution of disease is fixed by sampling). I also know that the odds ratio however can be calculated in such designs because it takes the same value when using the conditional distribution of $X|Y$ or $Y|X$. My question is: is it appropriate to use logistic regression in this case to model the odds of disease? i.e. $\text{logit}\left\{\dfrac{P(Y=1)}{1-P(Y=1)}\right\} = \mathbf{\beta} X $ context: GWAS (Genome Wide Association Studies) are typically case-control studies, where one wants to assess the association between disease and number of minor alleles of a particular SNP. $P$-values are typically obtained from a chi-squared test of independence. However, this doesn't allow you to add in other covariates. A lot of the packages that offer GWAS analysis also allow you to do logistic regression. I just wanted to verify that it was valid to do such an analysis.


title: Can I use an automated model selection approach on an lmer object?

text: I am attempting to use `MuMIn` to run a model selection analysis on a mixed model fitted using `lme4`. Because this model is fit using restricted maximum likelihood (REML), I am confused as to whether I can use an automated model selection approach such as that implemented in `MuMIn`. This is because `MuMIn` ranks models using AICc and maximum likelihood (ML), and thus the weights of any candidate models obtained are ML based. Does anyone know if this is a problem? If it is a problem, are there any other model selection packages in R that are more compatible?


title: Can Goldman Sachs' economists predict who will win the World Cup?

text: Goldman Sachs decided to apply their brains to predicting the winner of the World Soccer cup. They built a statistical model analyzing past winners. I read about it here, but the original version is in here. But I wonder why would results of past world cups be relevant at all? It is "a study of all mandatory international football matches since 1960 taking into account goals scored and conceded." I wonder whether including past results up to 1960 is just a joke or a blunder (similar to their other blunders). Otherwise, how can you analyze this statistically?


title: How to determine the conversion path?

text: A user interacts with X number of features before they go on to convert. How do I determine which features, **and in what order** , is the most common path to conversion? For example, if I was doing this analysis for Amazon I would have a "viewed product page" feature and an "add to wishlist" feature. My conversion would be a purchase. I would like an output that says something like: The most common path to conversion is when a user views 3 items, adds something to a wishlist, and then views 1 more item. I want to know the best order of events and the number of interactions with each feature. A decision tree, or what JMP calls Partition analysis, looks like it would be the perfect tool, but it doesn't take the order of events into account. Logistic regression is similar. It only takes the fact that an event happened into account and ignores the order in which the events occurred. How do I do an analysis that gives me insight into the order of inputs and not just the presence of inputs?


title: Conjugate mixture of beta distributions - weights

text: Suppose that my prior estimate for the bias of a coin is split evenly between a beta distribution (Theta,20,10) and a beta distribution (Theta,20,20), with Theta indicating the probability of tossing heads. Hence my prior is 0.5*Beta(20,10)+0.5*Beta(20,20). Assume that i observe that among 30 coin tosses i get 20 heads. I know that the updated betas have the form of beta(40,20) and beta(40,30), and that the new posterior has the form lambda*beta(40,20)+(1-lambda)*beta(40,30). however - i was hoping whether one can help me to understand how to calculate lambda / the weights of the updated distributions. Many thanks, W


title: Cross Validation (error generalization) after model selection

text: Note: Case is n>>p I am reading Elements of Statistical Learning and there are various mentions about the "right" way to do cross validation( e.g. page 60, page 245). Specifically, my question is how to evaluate the final model (without a separate test set) using k-fold CV or bootstrapping when there has been a model search? It seems that in most cases (ML algorithms without embedded feature selection) there will be   1. A feature selection step   2. A meta parameter selection step (e.g. the cost parameter in SVM). My Questions:   1. I have seen that the feature selection step can be done where feature selection is done on the whole training set and held aside. Then, using k-fold CV, the feature selection algorithm is used in each fold (getting different features possibly chosen each time) and the error averaged. Then, you would use the features chosen using all the data (that were set aside) to train the final mode, but use the error from the cross validation as an estimate of future performance of the model. _IS THIS CORRECT?_   2. When you are using cross validation to select model parameters, then how to estimate model performance afterwards? _IS IT THE SAME PROCESS AS #1 ABOVE OR SHOULD YOU USE NESTED CV LIKE SHOWN ON PAGE 54 (pdf) OR SOMETHING ELSE?_   3. When you are doing both steps (feature and parameter setting).....then what do you do? complex nested loops?   4. If you have a separate hold out sample, does the concern go away and you can use cross validation to select features and parameters (without worry since your performance estimate will come from a hold out set)? Many thanks to anyone who can clear this up! Sorry for a wordy question but I wanted to make sure I explained the confusion points (many ;-)) Brian


title: Analysis of variance table for 2x2 cross-over trials

text: I am trying to reproduce the anova table for a 2x2 cross-over design. I used the data listed in tables 2.1 and 2.2 of the "Design and Analysis of cross- over trials" book, by Jones and Kenward. I am giving all the code so that you can reproduce the calculations.               g1AB <- read.table(textConnection("     Label Per1 Per2     7 121.905 116.667     8 218.5 200.5     9 235 217.143     13 250 196.429     14 186.19 185.5     15 231.563 221.842     17 443.25 420.5     21 198.421 207.692     22 270.5 213.158     28 360.476 384     35 229.75 188.25     36 159.091 221.905     37 255.882 253.571     38 279.048 267.619     41 160.556 163     44 172.105 182.381     58 267 313     66 230.75 211.111     71 271.19 257.619     76 276.25 222.105     79 398.75 404     80 67.778 70.278     81 195 223.158     82 325 306.667     86 368.077 362.5     89 228.947 227.895     90 236.667 220     ")->con,header=T)     close(con)          g2BA <- read.table(textConnection("     Label Per1 Per2     3 138.333 138.571     10 225 256.25     11 392.857 381.429     16 190 233.333     18 191.429 228     23 226.19 267.143     24 201.905 193.5     26 134.286 128.947     27 238 248.5     29 159.5 140     30 232.75 276.563     32 172.308 170     33 266 305     39 171.333 186.333     43 194.737 191.429     47 200 222.619     51 146.667 183.81     52 208 241.667     55 208.75 218.81     59 271.429 225     68 143.81 188.5     70 104.444 135.238     74 145.238 152.857     77 215.385 240.476     78 306 288.333     83 160.526 150.476     84 353.81 369.048     85 293.889 308.095     99 371.190 404.762     ")->con,header=T)     close(con)          n1 <- nrow(g1AB)     n2 <- nrow(g2BA)     p <- 2 # periods      Some quantities that may be useful:               y11. <- sum(g1AB$Per1)     y12. <- sum(g1AB$Per2)     y21. <- sum(g2BA$Per1)     y22. <- sum(g2BA$Per2)     y1.. <- sum(c(g1AB$Per1,g1AB$Per2))     y2.. <- sum(c(g2BA$Per1,g2BA$Per2))     y... <- y1.. + y2..          y11.bar <- 1/n1 * y11.     y12.bar <- 1/n1 * y12.     y21.bar <- 1/n2 * y21.     y22.bar <- 1/n2 * y22.          y1..bar <- 1/(p*n1)*(y11.+y12.)     y2..bar <- 1/(p*n2)*(y21.+y22.)          y...bar <- 1/(p*(n1+n2)) * (y1.. + y2..)      In order to perform the analysis of variance in R, I created the following dataset               mydata1 <- data.frame(PEFR=g1AB$Per1,Subjects=g1AB$Label,Time=1,Groups="AB",Treatment=1)     mydata2 <- data.frame(PEFR=g2BA$Per1,Subjects=g2BA$Label,Time=1,Groups="BA",Treatment=2)     mydata3 <- data.frame(PEFR=g1AB$Per2,Subjects=g1AB$Label,Time=2,Groups="AB",Treatment=2)     mydata4 <- data.frame(PEFR=g2BA$Per2,Subjects=g2BA$Label,Time=2,Groups="BA",Treatment=1)     mydata <- rbind(mydata1,mydata2,mydata3,mydata4)     mydata$Subjects<-factor(mydata$Subjects)      (I have read that) The correct anova table is obtained using the following command               res <- summary(PEFR.aov <- aov(PEFR~Groups+Time+Treatment+Error(Subjects), data=mydata))          Error: Subjects               Df Sum Sq Mean Sq F value Pr(>F)     Groups     1  10573   10573   0.899  0.347     Residuals 54 634866   11757                         Error: Within               Df Sum Sq Mean Sq F value  Pr(>F)        Time       1    480   479.6   1.470 0.23061        Treatment  1   3026  3026.1   9.276 0.00359 **     Residuals 54  17617   326.2       1) Can you please explain to me the use of the Error term in the formula? 2) The anova table is the same to that of table 2.10 of the book, except the SS for Time (Period) which has a value of 396.858 instead of 480. According to Table 2.8 this quantity is calculated using               > n1*n2/(2*(n1+n2))*(y11.bar-y12.bar+y21.bar-y22.bar)^2     [1] 396.8583      Also, the Treatment SS is calculated by               n1*n2/(2*(n1+n2))*(y11.bar-y12.bar-y21.bar+y22.bar)^2     [1] 3026.12      If the formula was specified like this               PEFR~Groups+Treatment+Time+Error(Subjects)      the SS for Time would be correct, but the Treatment SS would have a value of 3109. How can I get 3026 for Treatment and 397 for Time, as I see in the book? What should be taken care of in the order the variables are specified in the formula? Thank you for your time.


title: Maximum likelihood estimator of SNR for a known signal superimposed in AWGN

text: Assuming the following model (Fir instant, received signal of a communication / RADAR system): $$ x[n] = \alpha s[n] + w[n] $$ Where:   * $ s[n] $ is a known signal of length $ L $ and known energy $ p = \sum_{n=0}^{L-1} {s}^{2}[n] $.   * $ \alpha $ is unknown attenuation factor. Though it is known $ \alpha \leq 1 $.   * $ w[n] $ is additive white gaussian noise (AWGN) where $ E[w[n]] = 0 $ and $ E[w[m]w[n] = \delta[n-m] {\sigma}^{2} $. The noise variance $ {\sigma}^{2} $ is unknown parameter.   * There are $ N $ samples, namely, $ n \in [0, \ldots, N-1] $ where $ N > L $. The question is, what is the Maximum Likelihood Estimator of the SNR of the received signal - $ \frac{{\alpha}^{2}p}{{\sigma}^{2}} $? Using Matched Filter / Correlation one could calculate the ML estimator of $ \alpha $ and $ {\sigma}^{2} $. Yet since the SNR isn't one to one function of the parameters, it cannot be calculated using them. Hence the question is: what's the MLE of the SNR of the received signal?


I found this book "Statistical methods in atmospheric sciences" by Daniel Wilks in the library. Do you think it is good for learning statistical methods in general? The examples are as you might know or expect, from atmospheric sciences but discussion about the statistical methods is quite good.


title: t-test and t-distribution etymology

text: Okay, I'm aware of the fact that Gosset came up with _t_ distribution, but what is the etymology of "t". How did "t" end up in t-test and t distribution?


title: Using the BH method for FDR control

text: I have preformed multiple tests and now have a list of all the calculated P-Values (~30K^2 values). I have a reason to believe my tests are somehow dependent. I would like to analyze the results and correct the p-values. Unfortunately, I only have very basic knowledge in statistics. From the little I managed to gather, I think the Benjamini Hochberg method for FDR control should be suitable for my needs. My questions are:   1. Am I right?   2. If so, I'd like to implement a code that preforms this calculation. Is the following code correct: http://www.augix.com/wiki/index.php/Benjamini_and_Hochberg_FDR_control and does it calculate the cutoff for my case? Thanks!!!


title: Major directions on a 2D map

text: Having a map such as:   ![enter image description here](http://i.stack.imgur.com/0NLWV.png)   **How to find major directions of variation** i.e., for the example given here: north-east--south-west? in an elegant, quick and clean way? Any ideas are more than welcome.


title: Advice on scientifically sound scale construction

text: I've been given a set of 20 Likert-items (ranging from 1-5, sample size n = 299) within the field of organizational research. The items are intended to measure a latent concept which is multidimensional, multifaceted and heterogenous in it's very nature. The goal is to create a scale(s) which can be nicely used to analyse different organizations and be used in logistic regression. Following the american association of psychology, a scale should be (1) unidimensional, (2) reliable and (3) valid. Therefore we decided to select four dimensions or subscales with 4/6/6/4 items each; which are hypothesized to represent the concept. The items werde constructed using the reflective approach (generating lots of possible items, and iteratively deleting items using cronbach's alpha and concept representation (validity) in three subsequent groups). Using the available data, a preliminary parallel explanatory factor analysis based on polychoric correlations and using varimax rotation revealed that items load on other factors than expected . There are at least 7 latent factors as opposed to four hypothesized ones. The average inter-item correlation is quite low (r=0.15) albeit positive. The cronbach-alpha coefficient is also very low (0.4-0.5) for each scale. I doubt that an Confirmatory factor analysis would yield a good model fit. If two dimensions were dropped, the cronbachs alpha would be acceptable (0.76,0.7 with 10 items per scale, which still could be made bigger by using the ordinal version of cronbachs alpha) but the scales themselves would still be multidimensional! As I'm new to statistics and lack the appropriate knowledge, I am at a loss on how to proceed further. As I am reluctant to discard the scale(s) completely and resign to a descriptive-only approach, I've got different questions: I) Is it wrong to use scales which are reliable, valid but not unidimensional? II) Would it be appropriate to interpret the concept afterwards as formative and use the vanishing tetrad test to assess model-specification and use partial least squares (PLS) to arrive at a possible solution? After all, the concept seems to be more a formative than a reflective one. III) Would using the item response models (Rasch, GRM etc) be of any use? As I've read, the rasch-models etc. also need the assumption unidimensionality IV) Would it be appropriate to use the 7 factors as new "subscales"? Just discard the old definition and use a new one based on factor loadings? I'd appreciate any thoughts on this one :) EDIT: Added factor loadings & correlations               > fa.res$fa     Factor Analysis using method =  ml     Call: fa.poly(x = fl.omit, nfactors = 7, rotate = "oblimin", fm = "ml")       Factor loadings calculated from factor pattern matrix and factor intercorrelation matrix, only values above 0.2 are displayed ![Factor Loadings](http://i.stack.imgur.com/9Ren4.jpg)


title: Which test is better suited to compare averaged versus single data sets

text: I have a light curve, that is photon count rate versus time, of an astronomical object. These data are periodic, since the source signal is periodic. I can fold the data with the period and obtain an average profile. At this point, I want to check whether the single profiles are statistically different from the averaged one. Then, I was thinking to do a test between the average profile and the single profiles. This test could be a 2sets KS test, or a Chi-squared test, or others. I don't know which one is better suited for my purpose. Could you suggest/address me in this please? EDIT: here is a figure where the folded profile is showed in blue points, and the red line a tentative of smoothing the profile (not too relevant for this discussion) ![enter image description here](http://i.stack.imgur.com/5PXgG.png) Also, here a small collection of "single-profiles" compared with the model from the red line in the previous figure http: fig1 , fig2 , fgi3. EDIT 2: here is a cumulative plot (very preliminary).


title: Missing data at random

text: How does one tell if a dataset is missing data at random? I've been reading up on how to impute missing values, and was wondering what techniques can be used to tell if data is really missing at random or systematically.


title: Interaction suppresses the main effect? How to interpret it?

text: I have a simple model without interaction and it stated significant effect for all the explanatory variables (continuous variable rok and categorical variables obdobi (levels hn and nehn) and kraj:               Call:     glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj, family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.8007  -1.1716  -0.5117   1.0864   4.2184            Coefficients:                   Estimate Std. Error t value Pr(>|t|)        (Intercept) -107.60761   53.96993  -1.994  0.04674 *      rok            0.05381    0.02686   2.003  0.04576 *      obdobinehn    -0.26962    0.11646  -2.315  0.02104 *      krajJHC        0.68869    0.31009   2.221  0.02683 *      krajJHM       -0.26607    0.32166  -0.827  0.40855        krajLBK       -1.11305    0.61942  -1.797  0.07298 .      krajMSK       -0.61390    0.41828  -1.468  0.14285        krajOLK       -0.49704    0.36981  -1.344  0.17958        krajPAK       -1.18444    0.39401  -3.006  0.00279 **     krajPLK       -1.28668    0.49672  -2.590  0.00988 **     krajSTC        0.01872    0.31222   0.060  0.95220        krajULKV      -0.41950    0.69220  -0.606  0.54478        krajVYS       -1.17290    0.44614  -2.629  0.00884 **     krajZLK       -0.38170    0.40969  -0.932  0.35198        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.645035)              Null deviance: 1136.22  on 489  degrees of freedom     Residual deviance:  938.02  on 476  degrees of freedom     AIC: NA          Number of Fisher Scoring iterations: 4      Then I added interaction obdobi:kraj:               Call:     glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj + obdobi:kraj,          family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.4635  -1.1706  -0.4597   1.0275   4.6829            Coefficients:                           Estimate Std. Error t value Pr(>|t|)        (Intercept)         -101.49501   54.53576  -1.861  0.06336 .      rok                    0.05102    0.02715   1.879  0.06086 .      obdobinehn            -1.11653    0.62058  -1.799  0.07264 .      krajJHC               -0.16805    0.51957  -0.323  0.74651        krajJHM               -0.77451    0.53738  -1.441  0.15018        krajLBK               -3.29567    1.42164  -2.318  0.02087 *      krajMSK               -0.73640    0.67267  -1.095  0.27420        krajOLK               -0.41582    0.68758  -0.605  0.54564        krajPAK               -1.50156    0.63871  -2.351  0.01914 *      krajPLK               -1.48611    0.75745  -1.962  0.05036 .      krajSTC               -0.34170    0.52059  -0.656  0.51191        krajULKV              -1.72550    1.02726  -1.680  0.09369 .      krajVYS               -1.93603    0.65862  -2.940  0.00345 **     krajZLK               -0.71065    0.65791  -1.080  0.28063        obdobinehn:krajJHC     1.44638    0.65507   2.208  0.02773 *      obdobinehn:krajJHM     0.82070    0.67910   1.209  0.22746        obdobinehn:krajLBK     3.31340    1.61026   2.058  0.04018 *      obdobinehn:krajMSK     0.12470    0.87281   0.143  0.88645        obdobinehn:krajOLK     0.04528    0.82529   0.055  0.95627        obdobinehn:krajPAK     0.48978    0.81921   0.598  0.55022        obdobinehn:krajPLK     0.23075    1.02316   0.226  0.82167        obdobinehn:krajSTC     0.50339    0.65976   0.763  0.44585        obdobinehn:krajULKV    2.49157    1.43679   1.734  0.08356 .      obdobinehn:krajVYS     1.48201    0.92082   1.609  0.10820        obdobinehn:krajZLK     0.49357    0.85087   0.580  0.56214        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.613648)              Null deviance: 1136.22  on 489  degrees of freedom     Residual deviance:  899.28  on 465  degrees of freedom     AIC: NA          Number of Fisher Scoring iterations: 4      **Strange thing happened - the main effects rok and obdobi are no longer significant! How can this happen? How to interpret this fact?** If the interaction `obdobi:kraj` has significant effect, then the `obdobi` also has significant effect, right? Note that the second model differs significantly (tested by `anova(..., test = "Chi")`). Thanks in advance! **EDIT:** added anova tables of the models (but since this is `glm` and not simple `lm`, mean sum of squares and p-values are missing and I don't know how to interpret it...)               > anova(model1)     Analysis of Deviance Table          Model: quasibinomial, link: logit          Response: cbind(ml, ad)          Terms added sequentially (first to last)                  Df Deviance Resid. Df Resid. Dev     NULL                      489    1136.22     rok      1     3.06       488    1133.16     obdobi   1    11.20       487    1121.96     kraj    11   183.94       476     938.02          > anova(model2)     Analysis of Deviance Table          Model: quasibinomial, link: logit          Response: cbind(ml, ad)          Terms added sequentially (first to last)                       Df Deviance Resid. Df Resid. Dev     NULL                           489    1136.22     rok           1     3.06       488    1133.16     obdobi        1    11.20       487    1121.96     kraj         11   183.94       476     938.02     obdobi:kraj  11    38.74       465     899.28


title: Response variable: percentage and too many zeros (zero inflated Poisson?)

text: I am analysing the effect of density (categorical), gonad mass (continuous) and temperature (continuous) on the percentage of acini spawning in a gonad. My replicate unit is a scallop. As my response variable is a percentage, and I have many zeros, I was wondering if a zero inflated Poisson regression would be adequate. I know this is used for count data - but as my data is percentage I am wondering if it would be all right to use this model.


title: A question about Latent Dirichlet Allocation model

text: when I used LDA model in my project, the result topic terms vary with the random seed. how to solve this problem ? thanks


title: A question on SEM and path analysis

text: For my PhD dissertation, I am offering a conceptual model on teacher professional development. My model has one mediator. I wonder if I can test the mediation through Baron and Kenny's procedure or bootstrap and then conduct structural equation modeling (SEM) analysis. Many say the first step (testing mediation) is an extra procedure, and they say SEM will handle mediation analysis. So, how can I justify using both procedures together? My second question is: should I propose a general research question for my PhD dissertation, such as, "Does $x$ mediate the relationship between exogenous and endogenous variables?" Or should I offer distinct research questions for every path in the model? I really appreciate your insightful responses.


title: Implementation of semi parametric methods

text: Has anyone worked with semi parametric methods to estimate parameters with binary outcome? Examples are like Cosslett (1983) or Ichimura or Klein-Spady. In other words we are looking for semi parametric methods instead of probit/logit. I was wondering if anyone has written R code or Matlab code in this regard that he/she could share?


title: Test which method is better if data is gathered on overlapping time windows (not independent)

text: I have two methods which are being used to estimate a specific signal. I have a ground truth measurement of this signal and these two methods are using noisy data to estimate this signal. This signal is computed using a sliding window, meaning that two consecutive estimates are not independent (the window is 30 seconds and a new estimate is made every 1 second, so the overlap in consecutive estimates is 29 seconds). So I now have two sets of estimates - one for each method. I want to determine which method is better, i.e. which one will have smaller residuals. I am using the absolute value of these residuals for the significance test, so the test has to be non-parametric (sort of a folded normal distribution). The two sets of estimates are paired - they are measured at the same point in time, just using different methods. I have been looking at Wilcoxon signed- rank test but there is an assumption that each pair is independent from each other, which they are not in this case. What is the best way to perform a statistical test on paired data for which each population is not independent? Please let me know if any clarification is required. Thanks in advance!


title: Visualising a linear model using effects package in R

text: I have ran this model in R:               model <- lm(mpg ~ hp + drat + disp, data=mtcars)      And I have visualised this model using `effects` package in `R`:               library(effects)     model_effects <- effect("hp", model, list(drat=mean(mtcars$drat), disp = mean(mtcars$disp)))     plot(model_effects)      ![enter image description here](http://i.stack.imgur.com/WHAZZ.jpg) I believe this plot shows the relationship between `mpg` and `hp` whilst holding `drat` and `disp` constant at their mean values. Could somebody confirm this is what the plot is doing?


title: Splines, latent variables and model identification

text: For $J$ units of observations my basic model contains latent variables $\theta_j, j = 1, ..., J$ identified by the marginal assumption $\theta \sim N(0,1)$. The latent variables are connected to a vector of observed binary variables $(x_{j1},...,x_{jI})$ via a binomial model and the logistic link: $logit(P(X_{ji} = 1)) = \alpha_i(\beta_i - \theta_j)$ (This is the 2PL or Birnbaum model, and the data is the matrix $(x_{ji})_{ji}$). Note that the latent parameters $\alpha_i$, $\beta_i$ are identified by the marginal assumption for $\theta$. I guess though my question could arise in other latent variable models. I now wish to extend the model, by including a metric covariate $y_j, j = 1,...,J$ that predicts $\theta_j$. This is not difficult to do for something of the form $\eta y_j$ or transformations of $y_j$. I tried using a spline for the influence of $y_j$ on $\theta_j$ and to ensure identifiability of the model, the spline should not include an intercept. My current approach to have an identified model is fix all members of the spline basis at a certain point, i.e. $R(\bar y, y_k^*) = 0$, $k = 1, ..., K$ here $\bar y$ is the mean of the $y_j$, $y_k^*$ is a knot of the spline, and $R(\cdot, z)$ is the spline basis function for the knot $z$. After reading parts of Wahba's (1990) and Woods (2006) book, I also realized that some spline bases satisfy constraints like $\int_0^1 R(y,z) dy = 0$, $\forall z$ which would lead to another approach w.r.t. identification. Curiously the first constraint (i.e. a spline with $f(\bar y) = 0$) seems to work, while the second seems unstable (i.e. unidentified). I have two questions now: 1) Are there any references on splines and latent variables that I might have missed? (I know there are people using splines for ICC estimation and estimation of marginal distributions since the 90s, but I could not find anything for covariates). 2) Is there a certain class of splines that is especially useful for this type of identification problem? (I have tried a truncated cubic basis and a basis given in Wood's book).


title: Calculating standard errors for contrasts of parameter estimates of a (multiple) linear regression

text: I am running regression analyses in matlab with glmfit and calculate significance of contrasts of parameter estimates by using standard errors. Somehow however, the results I get are not consistent with SPSS results, so I'm wondering if I am doing something wrong. Hence 2 Questions:   1. Can somebody point me to a source where the equations/formulas for of t-tests and F tests for contrasts of parameter estimates is described?   2. It would also be great if somebody could have a look of my code, to see if there is maybe something obviously wrong (see below). ### code Note: as this is matlab code, matrix algebra applies. 1) Calculate sigma square from residuals and degrees of freedom (stats.residp is a vector): `sigma_squared = (stats.residp'*stats.residp)/stats.dfe;` 2) Get the variance for a contrast through multiplication with covariance of regressors (stats.covb) and contrast vector (c, which is actually a matrix with many contrasts): `var_c = sigma_squared*c'*stats.covb*c;` 3) Get the standard error by taking the square root of the diagonal of the matrix of contrast-variances: `se_c = sqrt(diag(var_c));` 4) Calculate t-statistic (c,b,se_c are vectors; ./ makes the division an element by element division): `t = (c'*b)./se_c;` 5) Get the p value by plugging the t statistic in the cumulative student t distribution with the right degrees of freedom: `p = 2*((1-tcdf(t,stats.dfe)).*(t>0) + (tcdf(t,stats.dfe)).*(t<0));`


title: Meaning of hypothesis test for "μ = 25"; isn't it impossible?

text: I'm working through a stats textbook and have a question of the form: > You will perform a significance test of $H_0: μ=25$ based on an SRS of > $n=25$. Assume $σ=5$. I'm stuck on the 'equals' part of $μ=25$. Given that we're working in the set of reals, isn't the probability of $μ$ _equaling_ 25 zero? (Assume a non- textbook dataset; e.g. the odds of $μ=180$ for the height of adult humans is zero). I think what the hypothesis is trying to say is "$μ$ is within bounds dictated by your confidence level, but mean$(μ)=25$". The alternate hypothesis is something like $H_a: μ \neq 25$, which I know is true for any non-hypothetical dataset. So, what does $H_0: μ=25$ actually mean?


title: Statistical Test for Differences in Groups Based on Multivariate Function/Map

text: I've recently run into an issue where I am trying to test for a statistical difference between two groups, where each element of a group is itself a data object. For any pair of objects, I can calculate a similarity score. I want to check if the groups of objects are statistically different based on those similarity scores. Assume two groups of objects. Group A ($G_A^{\rightarrow}$) has $N_A$ objects and group B ($G_B^{\rightarrow}$) has $N_B$ objects. Let $x_{i,j}$ be object $j$ from group $i$, where $i$ is in {$A,B$} and $j$ is in {$0,1,...,N_i$} For any pair of objects, I have a function $f(x_{i_1,j_1},x_{i_2,j_2})$ that can give me a scalar score between 0 and 1 for the similarity of any two objects. So then, I can generate the following groups of scores:   * In-Group Similarities for Group A: $S_A$ = {$f(x_{A,j_1},x_{A,j_2}) \forall j_1,j_2 \in \{0,...,N_A,\}$}   * In-Group Similarities for Group B: $S_B$ = {$f(x_{B,j_1},x_{B,j_2}) \forall j_1,j_2 \in \{0,...,N_B,\}$}   * Between-Group Similarities: $S_{A * B}$ = {$f(x_{A,j_1},x_{B,j_2}) \forall j_1 \in \{0,...,N_A,\},j_2 \in \{0,...,N_B,\}$} I want to know if Group A has a significantly different distribution of objects, using similarity scores as a proxy. It is unclear what the appropriate method to approach this data would be. I either have two samples or three samples, depending on how you look at it. If I have two samples, then one is a combined sample of all in-group similarities (i.e., $S_A \cup S_B$) and the other is the between-group similarities (i.e., $S_{A*B}$). With two samples in hand, I could then hit them with some standard two-sample non- parametric tests. Alternatively, I could take each of these groups separately and end up with three samples and use tests for more than two groups (e.g., Kruskal-Wallis, Friedman, etc). In either case, it's a nightmare to interpret. It is easier to interpret a comparison of a single in-group sample against the between-group sample (e.g., $S_A$ vs $S_{A*B}$), but then similarity scores within the second group are completely ignored. Does anyone have experience applying traditional statistical tests to these? I'm trying to stick with them for interpretability. I could probably work out a Bayesian analysis, but that would be less ideal. Overall, there are three primary questions:   1. Can anyone recommend a good statistical test to apply for this problem?   2. Would a two sample (in-group vs between group), three sample (all separate), or other group structure be appropriate?   3. How would one interpret the probabilities generated by analyzing each of these samples? Notes:   * Obviously, I'd exclude duplicate pairings for calculating sets of similarity. I didn't bother representing this because it would be a pain.   * There are potentially multiple different similarity functions. You can't really rely on what the similarity function will look like, except that it is based on non-linear interactions between the objects (e.g., there is no concept of similarity without two objects).


title: Post hoc power analysis

text: I have an experiment with 3 equal group sizes and 4 measures. I think the simple null hypothesis is that the three groups will be the same. Most people, however, believe that group A should do best in measure 1, group B should do best in measure 2, group 3 should do best in measure C, and they should all be the same in measure 4, both when comparing within and between groups. It is quite possible my group size is too small to come to any conclusions. I base this on the tried and true approach that my N is smaller than their N. I started off by doing a multivariate GLM (using SPSS) and found that the effect group is not significant (F = 1.8, p = 0.1). Looking at the data suggests that even if there is a difference between the groups, that it doesn't likely follow what is expected (e.g., group A does not do the best in measure 1 and is best in measure 2). My thinking is that maybe a post hoc power analysis might tell me something useful. SPSS spits out an observed power of 0.69. At this point I am lost. Does the observed power tell me anything useful? Am I going about this all wrong?


title: Finding expectation of log

text: Does somebody have an idea to find the following expectation: $$E(X\log(a+bX)),$$ How can I proceed, can we go like this: $$=E(X)E[\log(a+bX)],$$if yes then what?


title: How do you perform a bayesian ordered logistic regression in R?

text: Trying to perform a Bayesian ordered logistic regression in R where age is my outcome variable. I have installed the ARM package but I am unsure how to go about generating my model in R. I also need to perform an ordered logistic in R using a Frequentist approach and use my model to generate predictive probabilities. Can you assist?


title: Variance-covariance matrix in lmer

text: I know that one of the advantages of mixed models is that they allow to specify variance-covariance matrix for the data (compound symmetry, autoregressive, unstructured, etc.) However, `lmer` function in R does not allow for easy specification of this matrix. Does anyone know what structure `lmer` uses by default and why there is no way to easily specify it?


title: Within model with plm package

text: I have a question about `plm` package. My code is:               fixedmodel <- plm(formula=Inv_ret~Firm.size+leverage+Risk+Liquidity+Equity,                       data=datanewp, model="within")      In the `plm vignette` the authors write: > This is called the fixed effects (a.k.a. within or least squares dummy > variables) model, usually estimated by ols on transformed data, and gives > consistent estimates for beta. I wonder then if `plm` package assumes that `least squares dummy variable` model is estimated with `model="within` in the `plm` package? In other words, if `wihin` is actually `least squares dummy variable` model in plm? On the other hand, if I refer to some other site they write: > There are several strategies for estimating a fixed effect model. The least > squares dummy variable model (LSDV) uses dummy variables, whereas the > “within” estimation does not. These strategies, of course, produce the > identical parameter estimates of regressors (nondummy independent > variables). So I guess the authors of `plm` package mean that they actually use `within` model but compute the unique group specific intercepts from `within` transformation later. Am I correct?


title: Multi-stage sampling together with hierarchical/ mixed effects models: which R packages?

text: Analyzing educational datasets we have samples of children from samples of class in samples of schools - we have sampling weights, so I use the survey package e.g. to do a linear model. But this kind of design also requires looking at the mixed effects. But this isn’t possible using the survey package. I can do this in nlme – but then I don't know how to account for the weighting. I guess I could use the sample weights as predictors in nlme regressions but I don’t think that is correct. It seems that this kind of design (in fact any stratified survey sample which includes nested levels) needs analysing from both perspectives – (survey weights and mixed effects) at once – but the packages of choice for each of these perspectives, survey and nlme, each don’t seem to have slots for the other perspective. Can someone put me on the right track, or suggest another package which does both at the same time?


title: Is there an algorithmic way to check if a dataset oscillates?

text: A simulation gives the population numbers for every species in the domain per frame. These vary over time and can be quite noisy, is there a way of analysing the data algorithmically without drawing a graph whether or not these values oscillate?


title: Choosing right algorithm for binary classification

text: I'm having a problem to predict hits from original features. I tried using LDA on original matrix but the thing is that probability of getting a hit vs non- hit is 95% vs 5%. That said after running LDA I get true positive rate of 5% and true negative rate of 96%. Now, my understanding of LDA makes me think that when it separates hits from non-hits into two groups the later gets more attention, however what I really want is to minimize the false discovery rate, even better to control it. The question is what algorithm is best suited for this? I have two binary coexclusive groups (Hit/Not-hit).


title: Figuring out probabilities with Hidden Markov Models

text: I'm really new to statistics so sorry in advance if this question does not make sense. **Background:** I'm trying to learn about hidden Markov models and they seem interesting but I was wondering about the probabilities they use to generate their predictions. Most of the information I have read has probabilities already about changing states and staying in the current state. I was on Github and looked up code(mainly python code) to understand a bit better but either they have the probabilities already computed or are calculating it based on occurrences in other documents(in the case of a spell checker, it studied a book to understand the words occurrences as well as the probability of them being beside each other). **Problem:** But I'm trying to understand the theory of figuring out probabilities for the HMM when I'm not using speech or text recognition. I'm trying to make a game which is trying to guess numbers based on a users behaviour(i.e. user is given a list of fruit prices and picks a basket of various fruits. We know the individual fruit prices and the users only reports the total cost of their basket, we then try to find the quantities of fruit they choose. They are not adversarial and cannot change their total quantities by more than 5%). My approach was to first find all the possibilities for all days, then using a graph to create edges between nodes that were within 5% and remove all improbable nodes. Now I have a range and with every turn new paths are being created and impossible ones are being removed. So within that range I am trying have a program 'zero' in on the correct solution. I think the HMM model would be helpful for this but I'm a bit confused on the application of Baum–Welch or forward-backward algorithm to create probabilities(that I can feed into my implementation of hmm). If it helps, I originally posted this at stackexchange (alot of details about the game itself as well code I have) . I got a very good solution, to find the nodes with the highest amount of edges as solutions. Its okay but it doesn't select a solution and if it makes a mistake it simply continues the path. I was also advised from mathexchange that this can be solved as a hidden markov model. I am basically struggling with applying it to my problem(right now figuring out the probability part. I have been told to create all paths as equally probable but I don't understand how that helps a HMM pick the most probable solution). If anyone has any suggestions or suggested reading or sample code(I'm learning python, but I can try to understand any other code thats provided). If possible can it be explained to a layman(I think I can apply the code if I can clearly understand the theory, which apparently I do not). Thanks in advance and sorry again if this question does not apply or is not explained correctly.


title: Understanding "variance" intuitively

text: What is the cleanest, easiest way to explain someone the concept of variance? What does it intuitively mean? If one is to explain this to their mom or child how would one go about it? It's a concept that I have difficulty in articulating - especially when relating variance to risk. I understand it mathematically and can explain it that way too. But when explaining real world phenomena how do you make one understand variance and it's applicability in the 'real world', so to speak. Let's say we are simulating an investment in a stock using random numbers (rolling a die or using an excel sheet, doesn't matter). We get some 'return on investment' by associating each instance of the random variable to 'some change' in the return. Eg.: > Rolling a 1 implies a change of 0.8 per \$1 in investment, a 5 a change of > 1.1 per \$1 and so on. Now if this simulation is run for about 50 times (or 20 or 100) we will get some values and the final value of the investment. So what does 'variance' actually tell us if we were to calculate it from the above data set? What does one "see" - If the variance turns out to be 1.7654 or 0.88765 or 5.2342 what does this even mean? What did/can I observe about this investment?? What conclusions can I draw - in lay man terms. Please feel free to augment the question with that for standard deviation too! Although I feel it's 'easier' to understand, but something that would contribute to making it also 'intuitively' clear would be greatly appreciated!


title: Sampling from kNN density estimate

text: Is there an easy/standard way of sampling from a density estimated by kNN? I've searched online but I can't seem to find how to do this. Any links or methods pointing me in the right direction would be much appreciated. Just to clarify, I would want to be able to sample points not (necessarily) in my dataset from the estimated density that I get using kNN on the dataset.


title: How is ARMA/ARIMA related to mixed effects modeling?

text: In time series analyses, I have used multi-level or random/mixed effects to deal with auto-correlation issues (i.e., observations are clustered within individuals over time) and added controls are added for some specification of time and for shocks of interest. ARMA/ARIMA seem designed to address similar issues. The resources I've found online discuss either ARMA/ARIMA or mixed effect models but I don't understand the relationship between the two. Might one want to use ARMA/ARIMA from within a multilevel model? Is there a sense in which the two are equivalent or redundant? Answers or pointers to resources that discuss this would be great.


title: Does using count data as independent variable violate any of GLM assumptions?

text: I would like to employ count data as covariates while fitting a logistic regression model. My question is:   * **Do I violate any assumption of the logistic (and, more in general, of the generalized linear) models by employing count, non-negative integer variables as independent variables?** I found a lot of references in the literature regarding hot to use count data as outcome, but not as covariates; see for example the very clear paper: "N E Breslow (1996) Generalized Linear Models: Checking Assumptions and Strengthening Conclusions, Congresso Nazionale Societa Italiana di Biometria, Cortona June 1995", available at http://biostat.georgiahealth.edu/~dryu/course/stat9110spring12/land16_ref.pdf. Loosely speaking, it seems that glm assumptions may be expressed as follows:   * iid residuals;   * the link function must correctly represent the relationship among dependent and independent variables;   * absence of outliers Does everybody know whether there exists any other assumption/technical problem that may suggest to use some other type of models for dealing with count covariates? Finally, please notice that my data contain relatively few samples (<100) and that count variables' ranges can vary within 3-4 order of magnitude (i.e. some variables has value in the range 0-10, while other variables may have values within 0-10000). A simple R example code follows:               \###########################################################          \#generating simulated data          var1 <- sample(0:10, 100, replace = TRUE);         var2 <- sample(0:1000, 100, replace = TRUE);         var3 <- sample(0:100000, 100, replace = TRUE);         outcome <- sample(0:1, 100, replace = TRUE);     dataset <- data.frame(outcome, var1, var2, var3);          \#fitting the model          model <- glm(outcome ~ ., family=binomial, data = dataset)          \#inspecting the model          print(model)          \###########################################################


title: Order of items in a questionnaire

text: Does the order of the items in a questionnaire matter? If so, how do I know which item should go first, second, third etc. when devising a questionnaire? Thanks


title: Techniques for scaling data matrix to avoid rank deficiency issues

text: I have a $n \times p$ matrix $A$ where $n$ is the number of observables and $p$ is the number of observations. $n \gg p$ In my code, I have done $[E,V] \,=\, eig(A)$ and doing a least squares operation on the eigenvector matrix $V$, which is $p \times p$. This is using the `mldivide` operator in MATLAB, such that `X = V\b;`. My data is fairly huge, and I am getting warnings from MATLAB that $V$ is rank deficient, and thus the condition number is too high to get a reasonable answer. _Are there any standard techniques/tools out there to overcome this?_ I have tried shifting and scaling $V$, and reading threads here on this topic, but no avail. The `mldivide` works for smaller datasets, where $p \sim O(10^1)$. However, when I go to $O(10^2)$ and beyond, I face rank deficiency issues. I would welcome suggestions. A few useful points about my data (I can provide more details if required):   1. This data is coming from an experiment of an extremely chaotic physical phenomenon (fluid turbulence), and therefore is noisy.   2. The $V$ matrix consists of eigenvectors in complex conjugates, and there is a lot of variation in every eigenvector (order of magnitude ranges from $10^{-14}$ to $10^{-1})$.


title: Estimating the asymptotic distribution of a quasi maximum likelihood estimator

text: We consider the following GARCH(1, 1) model: $y_t = h_t \epsilon_t$ where $(\epsilon_t)_{t \in \{1, \dots, n\}}$ are i.i.d. random variables with mean 0 and standard deviation 1. $h_t = \omega + \alpha \epsilon_{t-1}^2 + \beta h_{t-1}^2$ A standard way of estimating the parameters of such a model is to use the quasi maximum likelihood method, i.e. maximising the log-likelihood of the model under the assumption that the $(\epsilon_t)_{t \in \{1, \dots, n\}}$ are gaussian, that is: $\displaystyle\mathscr{l}(\omega, \alpha, \beta, y_1, \dots, y_n) = -\frac{1}{2}\log(2\pi) - \frac{1}{2}\sum_{t=1}^n\left(\log(h_t) + \frac{y_t^2}{h_t}\right)$ The estimator is shown to have good properties even if this assumption does no hold. The asymptotic distribution of that estimator is: $\sqrt{n}\left(\hat{\theta}_n - \theta_0\right) \rightarrow \mathcal{N}\left(0, A(\theta_0)^{-1} B(\theta_0) A(\theta_0)^{-1}\right)$ where $\hat{\theta}_n$ is the vector of estimated parameters, $\theta_0$ the true vector of parameters, $B(\theta_0)$ the outer product of the gradient and $A(\theta_0)$ the Hessian matrix, both evaluated at the true value of the parameters. We can get a good approximation of the Hessian asymptotically by replacing the true value of the parameter by its estimate. But when it comes to the gradient, $B(\hat{\theta}_n)$ is by definition equal to zero since $\hat{\theta}_n$ maximises $\mathscr{l}$. So it cannot be a good approximation of $B(\theta_0)$. That leads me to my question: how do you estimate the covariance matrix of the asymptotic distribution of the parameters?


title: A one-way between subjects ANOVA

text: I have 3 visualization techniques (e.g., PC, NL and Map). I wanted to evaluate user performance (Time and Accuracy) for these three visualizations methods. I have 9 tasks for participants to perform. These 9 tasks categorized into 3 groups such group 1 , group 2 and Group 3 (each groups contains 3 tasks). I ran between subjects with 36 participants (12 for each visualization method). I asked participants to perform same tasks in all three visualizations. Now I have 6 table which task completion time and Accuracy (whether participant could successfully perform the tasks or not) for 3 different visualization.               One:  Completion time for PC visualization method                     Participant 1 ......... Participant 12       Task 1     23s (task completion time)                    Task 9               Second : Accuracy for  PC visualization method                     Participant 1 ......... Participant 12       Task 1     Yes (yes if task completed)                         Task 9      How should I analyse these 6 tables? should I apply Anova? if yes which Anova method is the best one for my case? How should I apply Anova for accuracy where the results for accuracy is binary?


title: Statistical videos to learn regression analysis and multivariate analysis?

text: I'm having a tough time to understand concepts of Regression Analysis and Multivariate Analysis. I'm following the books of Johnston and Anderson, but some video lectures would help me learn it faster. Could anyone suggest me some videos of the same? Thanks.


title: How to express an hypothesis that something has increased?

text: I am trying to express a hypothesis to test for the following claim: > A biologist is presented with the data that shows an increase in the average > number of bacterias, though he suspects there was no actual change. I came up with the following:   * $H_0:\; \mu = 40.10^2$ [no change had occurred]   * $H_A:\; \mu > 40.10^2$ [there is an increase]


title: Probability of drawing no red balls from 20 draws without replacement given finite sample

text: I understand this to be a binomial distribution: There are 100 balls in a bucket. 10 are red, 90 are blue. I select a ball at random and then replace it in the bucket, and I do this 20 times. I then calculate the probability that none of the selected balls were red. But what if I don't put the ball back in the bucket? The probability then changes with each trial. Can anyone point me in the right direction of how to calculate the probability in this case?


title: Recursive partitioning using rpart() method in R

text: I am new to R and using `rpart` for building a regression tree for my data.I wanted to use all the input variables for building the tree, but the rpart method using only a couple of inputs as shown below. As we can see, I have provided 10 inputs, but rpart used only two inputs. Please let me know how can force rpart method to use all the input variables. Thanks.               rm = rpart(uloss ~ tc_b + ublkb + mpa_a + mpa_b +           sys_a + sys_b + usr_a, data = data81, method="anova")     > princtp(rm)            Regression tree:     rpart(formula = uloss ~ tc_b + ublkb + mpa_a + mpa_b + sys_a +          sys_b, data = data81, weights = usr_a, method = "anova")          Variables actually used in tree construction:     [1] mpa_a tc_b          Root node error: 647924/81 = 7999     n= 81                CP nsplit rel error  xerror     xstd     1 0.403169      0   1.00000 1.04470 0.025262     2 0.092390      1   0.59683 0.66102 0.015238     3 0.081084      2   0.50444 0.70702 0.013123     4 0.045304      3   0.42336 0.58683 0.012129     5 0.010000      4   0.37805 0.51930 0.011942      One more question: I have used rpart.control for minsplit=2, and got the following for another data. Inorder to avoid overfititng the data, do I need to use splits 3 or splits 7. Shouldn't I use splits 7? Please let me know. Variables actually used in tree construction: [1] ct_a ct_b usr_a Root node error: 23205/60 = 386.75 n= 60                       CP nsplit rel error  xerror     xstd     1 0.615208      0  1.000000 1.05013 0.189409     2 0.181446      1  0.384792 0.54650 0.084423     3 0.044878      2  0.203346 0.31439 0.063681     4 0.027653      3  0.158468 0.27281 0.060605     5 0.025035      4  0.130815 0.30120 0.058992     6 0.022685      5  0.105780 0.29649 0.059138     7 0.013603      6  0.083095 0.21761 0.045295     8 0.010607      7  0.069492 0.21076 0.042196     9 0.010000      8  0.058885 0.21076 0.042196


title: Problem with maximum likelihood estimation in R: NaNs produced

text: I'm trying to estimate 4 maximum likelihood estimators on a data set containing two variables: `x` and `y`. When using the following code:                   logLikGrad <- function(param) {         b0 <- param[1]         b1 <- param[2]         sigma <- param[3]         landa <- param[4]         logLikGradValues <- numeric(2)         logLikGradValues[1] <- sum( log(1+(y-b0-b1*x)/sigma*(x^landa)*3))         return(logLikGradValues)         }          mleGrad <- maxLik(logLik = logLikFun, grad = logLikGrad, start = c(b0 = 1, b1=2,sigma = 1, landa=1))      Using the code above gives me the following warning message: `"In log(1 + (y - b0 - b1 * x)/sigma * (x^landa) * 3) : NaNs produced"` Does anyone know how to solve this problem?


title: How to find groupings (trajectories) among longitudinal data?

text: ### Context I want to set the scene before somewhat expanding on the question. I have longitudinal data, measurements taken on subjects approximately every 3 months, primary outcome is numeric (as in continuous to 1dp) in the range 5 to 14 with the bulk (of all data points) being between 7 and 10. If I do a spaghetti plot (with age on the x axis and a line for each person) it's a mess obviously as I have >1500 subjects, but there is a clear tread towards higher values with increased age (and this is known). **The broader question:** What we would like to do is to firstly be able to identify trending groups (those that start high and stay high, those that start low and stay low, those that start low and increase to high etc) and then we can look at individual factors that are associated with 'trend group' membership. My question here is specifically pertaining to the first portion, the grouping by trend. ### Question   * How can we group individual longitudinal trajectories?   * What software would be suitable for implementing this? I have looked at Proc Traj in SAS and M-Plus suggested by a colleague, which I'm looking into, but would like to know what others thoughts are on this.


title: Help evaluating a posterior probability expression

text: Consider $\boldsymbol{x}= [x_1,x_2,...x_n]$ and $\boldsymbol{y}= [y_1,y_2,...y_n]$ to be two multivariate Gaussians with an isotropic diagonal variance structure and uninformative priors so that: $p(x_i,\mu_{x_i},\sigma_x) = \frac{e^{\frac{-(x_i-\mu_{x_i})^2}{2\sigma_x^2}}}{\sqrt{2\pi \sigma_x}}$ and, $p(y_i,\mu_{y_i},\sigma_y) = \frac{e^{\frac{-(y_i-\mu_{y_i})^2}{2\sigma_y^2}}}{\sqrt{2\pi \sigma_y}}$ where $\sigma_x$ is a scalar, constant for all $i$, $\sigma_y$ is a scalar, constant for all $i$; and $\mu_{x_i}$ and $\mu_{y_i}$ are scalars for each $i$. We can define vectors $\boldsymbol{\mu_x} = [\mu_{x_1},...\mu_{x_n}]$, and $\boldsymbol{\mu_y} = [\mu_{y_1},...\mu_{y_n}]$. For simplicity, let us assume that we have unnormalized uninformative uniform priors of all $\mu$ and $\sigma$ so that, $p(\boldsymbol{x},\boldsymbol{\mu_x},\sigma_x) =\Pi_i p(x_i,\mu_{x_i},\sigma_x) = \Pi_i \frac{e^{\frac{-(x_i-\mu_{x_i})^2}{2\sigma_x^2}}}{\sqrt{2\pi \sigma_x}}$ and $p(\boldsymbol{y},\boldsymbol{\mu_y},\sigma_y) =\Pi_i p(x_i,\mu_{x_i},\sigma_x) = \Pi_i \frac{e^{\frac{-(y_i-\mu_{y_i})^2}{2\sigma_y^2}}}{\sqrt{2\pi \sigma_y}}$ Given this, can I evaluate $p(\boldsymbol{\mu_y} \propto \boldsymbol{\mu_x}|\boldsymbol{x},\boldsymbol{y})$ i.e., $p(\boldsymbol{\mu_y}=m\boldsymbol{\mu_x}+c|\boldsymbol{x},\boldsymbol{y})$ for any scalars $m$ and $c$? Just to explain what I mean, since we already have expressions for $p(\boldsymbol{x},\boldsymbol{\mu_x},\sigma_x)$ and $p(\boldsymbol{y},\boldsymbol{\mu_y},\sigma_y)$, I am interested in $\frac{\int p(\boldsymbol{x},\boldsymbol{\mu_x},\sigma_x)p(\boldsymbol{y},\boldsymbol{\mu_y},\sigma_y)\mathbf{1}_{[\frac{\mu_{y_1}-\bar{\mu_y}}{\mu_{x_1}-\bar{\mu_x}}=\frac{\mu_{y_2}-\bar{\mu_y}}{\mu_{x_2}-\bar{\mu_x}}=...=\frac{\mu_{y_n}-\bar{\mu_y}}{\mu_{x_n}-\bar{\mu_x}}]}d\sigma_x d\sigma_y d\boldsymbol{\mu_x} d\boldsymbol{\mu_y}} {\int p(\boldsymbol{x},\boldsymbol{\mu_x},\sigma_x)p(\boldsymbol{y},\boldsymbol{\mu_y},\sigma_y)d\sigma_x d\sigma_y d\boldsymbol{\mu_x} d\boldsymbol{\mu_y}}$ where $\bar{\mu_x} = \Sigma_i \mu_{x_i}/n$ and $\bar{\mu_y} = \Sigma_i \mu_{y_i}/n$, and $\mathbf{1}_{[P]}$ is the indicator function, $1$ if proposition $P$ is true, $0$ otherwise. I thought of doing this by removing $\mu_y$ in the numerator, replacing it by $m\mu_x+c$, and integrating over $m$ and $c$, but I do not know what pdf to use for $m$ and $c$ for the uninformative $\mu$s, surely $m$ and $c$ are not simply uniform.


title: Test of normality repeated measure ANOVA

text: I'm trying to do a repeated measure ANOVA analysis on the 5 NEO-FFI factors (Neuroticism, Extraversion, Openness to experience, Agreeableness and Conscientiousness). Sujest are tested only once, so I just want to compare the factors. My question is how do I do the normality test? I mean, what is the dependent variable?


title: Binary classification when many binary features are missing

text: I'm working on a binary classification problem, with about 1000 binary features in total. The problem is that for each datapoint, I only know the values of a small subset of the features (around 10-50), and the features in this subset are pretty much random. What's a good way to deal with the problem of the missing features? Is there a particular classification algorithm that handles missing features well? (Naive Bayes should work, but is there anything else?) I'm guessing I don't want to do some kind of variable imputation, since I have so many missing features.


title: One Factor Repeated Measure, Two Way ANOVA

text: Is there some way in R to output not only the main effects and the interactions of the two factors using an aov() like function but also the variance of the subject within a non-repeated factor? For example, I have an experiment dealing with the factor Dose(the between subject variable - non repeated) and Time(the within subject - repeated) on 40 animal's AUC in a balanced design (10 animals in each of 4 doses measuring AUC at 2 times, visit 1 and visit 2, for a total of 80 observations). The aov() function outputs the Dose, Visit, and Dose:Visit interaction just fine, but the residual output contains pertinent information I'd like to have. This includes the Animal w/i Dose Sums of Square as well as the within animal error, both of which have a df of 36. How would I obtain these statistics in R? Any help is appreciated, thanks.


title: How to handle cut-off weeks in time series analysis?

text: I am doing a time series analysis of sales data from a data warehouse. For that I want to use data grouped by the week of year. My problem now is that e.g. for week 1/2014, I have an outlier as this week only has two work days here in Austria. Is such a thing a problem? If yes, what is considered as good practice to handle this?


title: Optimal stopping under partially observable state

text: This problem is basically the classic asset selling problem but with imperfect state information. In the classical problem, we have an asset that we wish to sell, we receive offers w(0) to w(N-1). If we accept the offer at a given period, we can invest that money with a certain rate of interest r > 0. We suppose that at period N-1, we will accept the offer if we haven't accepted any already. In the problem I am trying to solve, the offers w(k) are independent and identically distributed. However the common distribution of the w(k) is unknown. Instead it is known that this distribution is one out of two given distributions F1 and F2, and that the a priori probability that F1 is the correct distribution is a given scalar q, with 0 < q < 1. a)Formulate this as an imperfect state information problem and identify the state, control, system disturbance, observation, and observation disturbance. b)Show that (x(k), q(k)), where q(k) = P(distribution is F1 | w(0), ... , w(k-1)) is a suitable sufficient statistic, write a corresponding DP algorithm, and derive the form of the optimal selling policy. The state would represent whether or not we have already sold our asset and if we haven't, the offer we have at this period. The control would be to either accept or refuse the offer at this period. The system disturbance would be the offer. As for the rest of the question, I'm not too sure.


title: Determining the nature of noise

text: I have a general question. What kind of noise is additive, and what about multiplicative noise? How to determine the nature of noise? Thanks a lot for your help.


title: Binomial testing with probability estimated from sample

text: If I have a sample of k successes and n-k failures, there are standard techniques (Agresti-Coull, Clopper, etc.) for finding a confidence interval of the probability of an individual success. What if I want to find a confidence interval for the probability of getting at least k' out of n' instead? Obviously it can be approximated by the usual $$\sum_{j\ge k'}{n'\choose j}(k/n)^j(1-k/n)^{n'-j}$$ but this does not take into account the uncertainty of the initial sample, which may be important (the samples are small).